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ABSTRACT 


The  U.S.  Army’s  Future  Force  requires  information  dominance  to  succeed,  yet 
finds  itself  with  an  ever-increasing  gap  between  its  capacity  to  collect  infonnation  and  its 
information  processing  capacity — with  little  understanding  of  how  to  efficiently  utilize 
scarce  processing  resources.  In  this  investigation,  a  model  is  proposed  to  adequately 
represent  the  flow  of  information  within  a  command  and  control  context  toward  the  end 
of  optimally  controlling  this  flow.  The  model  is  conjectured  to  be  NP-hard  in  general. 
Closed-fonn  optimal  solutions  are  derived  for  special  cases  of  the  model,  while  other 
cases  are  shown  to  be  NP-hard.  A  case  of  the  model  is  shown  to  equate  to  a  special  case 
of  the  quadratic  assignment  problem  not  previously  known  to  have  a  closed-form 
solution,  and  such  a  solution  is  derived.  Upper  and  lower  bounds  are  derived  for  more 
general  cases  of  the  model,  and  heuristic  strategies  are  proposed  and  tested  in  discrete 
event  simulation.  Strong  empirical  evidence  is  produced  to  demonstrate  the  effectiveness 
and  robustness  of  one  heuristic. 
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I.  INTRODUCTION 


A.  BACKGROUND  AND  MOTIVATION 

The  U.S.  Army’s  Future  Force  is  envisioned  to  “see  first,  understand  first,  act 
first,  and  finish  decisively”  [1].  “Seeing  first”  is  accomplished,  in  part,  through  the 
deployment  of  multiple  manned  and  unmanned  systems  throughout  the  operational 
environment.  Future  Army  commanders  are  expected  to  have  access  to  infonnation 
collected  from  an  enonnous  array  of  Intelligence,  Surveillance,  and  Reconnaissance 
(ISR)  assets.  Brigade  level  commanders — such  as  those  commanding  the  Army’s  Future 
Combat  System  Brigade  Combat  Teams  (FCS  BCT) — will  directly  control  hundreds  of 
individual  manned  and  unmanned  systems  and  will  have  access  to  infonnation  collected 
by  ISR  assets  belonging  to  sister  organizations  and  higher  echelons,  to  include  Joint  and 
strategic  assets.  The  FCS  BCT  is  programmed  to  have  a  total  of  8  different  ISR  systems 
(air,  ground,  manned,  and  unmanned)  in  addition  to  the  planned  combat  systems — which 
come  with  their  own  onboard  ISR  capabilities  [2],  Together,  these  assets  will  accomplish 
a  variety  of  collection  missions  in  support  of  the  commander’s  critical  information 
requirements  (CCIR)  and  intelligence  collection  plan,  to  include  wide  area  surveillance 
(WAS)  and  reconnaissance,  surveillance  and  target  acquisition  (RSTA).  In  addition  to 
the  ISR  assets  currently  planned  for  the  FCS  BCT,  new  technologies  will  undoubtedly  be 
introduced  that  will  increase  even  further  the  amount  of  information  available  to  the 
commander.  Some  of  these  technologies,  to  include  autonomous  UAV  swarms  [3]  and 
self-forming  networks  of  tiny  “smart  dust”  sensors  [4],  are  presently  emerging.  The  data 
and  infonnation  made  available  to  the  commander  combine  to  form  an  immense 
collection  of  information  packages  taking  many  possible  forms,  including  voice 
transmissions,  text  messages,  video,  still  imagery,  and  radar  tracking  data,  as  well  as  data 
from  magnetic,  acoustic,  seismic,  chemical/biological,  and  weather  sensors.  The  “see 
first”  component  of  the  Army’s  vision  appears  well  on  its  way  to  being  realized. 
“Understanding  first,”  however,  specifically  the  Anny’s  capacity  (or  relative  lack 
thereof)  to  process  and  prioritize  the  vast  amounts  of  intelligence  collected,  will  remain 
an  enonnous  challenge  into  the  foreseeable  future. 
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The  problem  of  finding  a  small  subset  of  relevant  information  in  a  very  large, 
rapidly  changing  data  set  is  not  new.  Data  mining  and  sensor  data  fusion  approaches 
have  been  used  recently  to  address  this  issue  (see  [5]  for  example).  This  research, 
however,  focuses  on  a  separate — yet  related — problem:  total  information  flow  in  support 
of  decision  making.  This  research  emphasizes  the  dynamic  nature  of  information  flow 
and  the  interaction  between  infonnation  processing  and  the  commander’s  decision 
making  process.  Here  we  introduce  a  dynamical  battlefield  information  flow  model.  The 
model  is  designed  so  that  it  is  able  to  accommodate  a  variety  of  important  factors  such  as 
information  volume  and  value,  efficiency  of  information  handling,  commander’s 
feedback,  and  random  factors.  We  first  introduce  an  existing  descriptive  model,  the 
Dynamic  Model  of  Situated  Cognition,  as  a  lead-in  to  the  dynamic  information  flow 
modeling  that  follows. 

The  Dynamic  Model  of  Situated  Cognition  (DMSC)  [6]  was  first  introduced  in 
2003  and  has  quickly  gained  acceptance  in  the  Department  of  Defense  community  as  a 
descriptive  model  of  information  processing  and  information  flow  on  the  battlefield. 
Shown  in  Figure  1,  this  model  illustrates  the  relationship  between  technological  systems 
and  human  perceptual  and  cognitive  processes  and  provides  an  effective  framework 
within  which  to  view  infonnation  flow  in  a  complex  system.  While  not  an  analytical 
model,  the  DMSC  does  provide  a  means  of  viewing  the  information  flow  problem  in  a 
way  that  is  accessible  to  analysts  and  decision  makers  alike. 

The  technological  portion  (the  left  side)  of  the  model  characterizes  information  as 
being  present  in  one  or  more  successively  filtered  “ovals.”  Oval  1  is  Ground  Truth, 
consisting  of  all  available  information  at  a  given  point  in  time.  Oval  1  is  a  completely 
accurate  description  of  ground  truth  and  is  constantly  updated  (and  therefore  dynamic). 
Oval  1  contains  infonnation  about  friendly  and  enemy  forces,  to  include  locations  and 
sizes  of  forces  and  available  weapons  systems,  as  well  as  less  tangible  infonnation  such 
as  commanders’  intent  and  troop  morale.  The  shapes  shown  in  Oval  1  represent 
individual  information  “packages”  pertaining  to  friendly  and  enemy  entities,  their 
disposition,  and  commanders’  intent  as  well  as  data  on  noncombatants,  weather,  terrain 
and  other  environmental  conditions. 


2 


All  data  in  the 
environment 


Projection  of 
decision  maker 

©  Miller  and  Shattuck,  2003 


Figure  1  Dynamic  Model  of  Situated  Cognition 

(from  [6]) 


Oval  2  is,  in  general,  an  error-prone  subset  of  the  information  contained  in  Oval  1 
and  represents  the  complete  set  of  information  that  has  been  detected  or  collected  by 
technological  systems  (i.e.,  sensors).  Two  types  of  error  are  found  in  Oval  2 — 
information  can  be  missed  entirely,  or  it  can  be  misrepresented.  Information  can 
generally  be  misrepresented  in  two  ways:  an  entity  can  be  represented  as  a  threat  when  it 
is  not  (e.g.,  a  decoy  or  a  noncombatant),  or,  conversely,  a  true  threat  can  be  represented 
as  a  non-threatening  entity.  Accurate  and  inaccurate  representations  are  characterized  by 
the  various  shapes  contained  in  Oval  2  in  Figure  1.  As  with  Oval  1,  Oval  2  is  dynamic — 
the  infonnation  contained  in  it  changes  with  time. 

Oval  3  contains  all  information  shown  on  a  decision  maker’s  command  and 
control  (C2)  display.  It  contains  a  filtered,  fused  subset  of  the  information  contained  in 
Oval  2.  The  future  commander  is  expected  to  have  the  ability  to  tailor  his  own  display  by 
accepting  various  levels  of  autonomous  configuration  in  combination  with  manual 
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configuration  controls  (zoom,  filtering,  field  of  view,  etc.).  Clearly,  any  error  present  in 
Oval  2  can  propagate  to  Oval  3  (and  beyond).  Furthermore,  additional  error  is  introduced 
through  inaccurate  fusion  algorithms  and  improper  filtering  techniques.  Again,  the 
representation  of  information  in  Oval  3  is  dynamic. 

The  three  remaining  ovals  comprise  the  cognitive  portion  of  the  model.  Oval  4 
contains  all  information  perceived  by  the  decision  maker.  Information  contained  in  Oval 
4  is  a  function  of  that  in  Oval  3,  but  not  (usually)  a  one-to-one  correspondence.  Oval  5 
represents  the  comprehension  of  the  decision  maker  based  on  perceived  infonnation  and 
his  own  individual  lens.  Finally,  Oval  6  represents  the  projection  or  prediction  of  the 
decision  maker. 

Preceding  each  of  the  last  three  ovals  is  a  lens  signifying  individual  filters  through 
which  information  passes  as  a  decision  maker  moves  from  seeing  to  perceiving,  then  to 
comprehending  and  finally  projecting.  These  lenses  are  based  on  the  individual’s 
experience  level,  training,  operational  orders  (OPORDs),  guidance  and  policy,  and 
localized  factors  such  as  fatigue  or  fear.  Lenses  can  be  skewed — resulting  in  errors  in 
perceiving,  comprehending,  and/or  predicting — based  on  errors  that  have  propagated 
from  Ovals  1-3  and  lack  of  experience  and/or  training,  or  fear,  fatigue  and  other  factors 
[6]. 

The  dynamic  nature  of  the  DMSC  is  seen  in  Figure  2.  Infonnation  can  be  thought 
of  as  flowing  from  one  oval  to  the  next,  and  decisions  are  reflected  as  feedback  within  the 
system,  dynamically  altering  the  contents  of  each  oval.  Feedback  can  take  the  fonn  of 
actions  in  the  real  world  (Oval  1),  redistribution  of  sensors  (Oval  2),  and  information 
filtering  on  local  C2  systems  (Oval  3).  For  example,  based  on  what  he  comprehends 
about  the  enemy  situation,  disposition,  and  intent,  the  commander  may  decide  to  engage 
the  threat.  This  will  cause  a  change  in  ground  truth  (Oval  1) — due,  perhaps,  to  enemy  or 
friendly  casualties,  or  enemy  displacement — and  cascading  changes  in  all  subsequent 
ovals  [7].  Another,  less  obvious,  example  of  feedback  in  the  system  takes  place  when  the 
commander  modifies  existing  information  processing  priorities.  This  will  result  in  a 
change  to  the  information  he  sees  on  his  C2  display  (Oval  3),  which  in  turn  affects  his 

perception,  comprehension  and  projection  (Ovals  4-6). 
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Figure  2  Representation  of  feedback  in  the  Dynamic  Model  of  Situated  Cognition 

(from  [8]) 

To  this  point,  we  have  described  a  system  within  which  information  on  the 
battlefield  dynamically  flows  and  changes  based  on  controllable  factors  (friendly 
decisions  and/or  actions)  as  well  as  uncontrollable  factors  (enemy  decisions/actions, 
random  factors).  In  the  following  section,  we  describe  a  family  of  mathematical  models 
that  demonstrate  the  potential  to  1)  accurately  portray  this  information  flow  system,  and 
2)  inform  the  process  by  which  information  flow  is  controlled  and  prioritized. 

B.  DEFINING  AND  SCOPING  THE  PROBLEM 

Figure  3  illustrates  several  potential  mathematical  models  and  the  portions  of  the 
DMSC  they  align  with.  First,  we  have  a  dynamic  sensor  allocation/reallocation  model 
(corresponding  to  Ovals  1  and  2  of  the  DMSC).  This  model  could  potentially  produce  an 
allocation  or  reallocation  plan  for  organic  sensor  platforms  based  on  the  following  inputs, 
among  others. 
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•  Information  on  organic  sensor  platforms — location,  disposition, 
capabilities,  range,  maximum  speed,  cruising  speed,  etc. 

•  CCIR — information  elements  that  directly  contribute  to  successful  mission 
accomplishment  [9] 

•  Likely  threat  disposition  and  intent 

•  Friendly  mission 


Figure  3  Family  of  Dynamical  System  Models 

Data  and  information  fusion  algorithms  and  models  can  also  be  related  to  the 
DMSC.  They  reside  in  the  space  between  Ovals  2  and  3,  taking  raw,  unprocessed  data 
and  producing  processed,  fused  and  fdtered  information  for  display  on  a  C2  system. 
Additionally,  human  cognition  models  (Ovals  3-6)  representing  naturalistic  decision 
making  and  other  approaches  would  be  appropriate  within  the  context  of  the  DMSC,  as 
would  models  that  represent  the  dynamics  of  a  commander’s  situational  awareness 
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(Ovals  1-5).  These  four  modeling  areas  (sensor  allocation,  fusion,  human  cognition,  and 
situational  awareness)  are  well-developed  in  the  literature  (see  [5,  10-13]  for  examples) 
and  are  mentioned  for  contextual  purposes  only.  They  are  beyond  the  scope  of  this 
research. 

This  research,  then,  focuses  on  the  modeling,  analysis,  and  simulation  of  dynamic 
information  flow  within  the  command  and  control  system  (the  physical  system)  as  shown 
in  Figure  3.  There  are  two  references  to  this  modeling  area  in  Figure  3:  modeling  the 
dynamics  of  information  flow  within  the  entire  system  remains  a  long-term  research  goal; 
the  true  focus  of  this  research,  however,  is  on  modeling  the  dynamical  flow  of 
information  between  Ovals  2  and  3.  Herein  lies  the  bridge  between  a  commander’s 
ability  to  see  first  and  understand  first.  This  modeling  area  will  be  developed  in  detail  in 
subsequent  chapters. 

C.  REVIEW  OF  PREVIOUS  RESEARCH 

Two  main  areas  in  the  literature  are  touched  upon  in  this  research:  the  family  of 
Quadratic  Assignment  Problems  (QAPs)  and  the  area  of  sequencing  and  scheduling.  We 
shall  address  each  of  these  in  this  literature  review. 

1.  Quadratic  Assignment  Problem 

Koopmans  and  Beckman  [14]  first  introduced  the  Quadratic  Assignment  Problem 
(QAP)  in  1957  as  a  mathematical  model  for  the  assignment  of  n  “indivisible  economic 
activities”  (i.e.,  plants)  to  n  locations.  The  general  QAP  is  known  to  be  NP-complete  [15, 
16].  Burkard,  et  al.  [17]  identify  special  cases  of  the  QAP  having  polynomial-time 
solutions.  Additionally,  Burkard,  et  al.  [17-20],  Cela  [21,  22],  and  most  recently 
Demidenko  [23]  in  2006  prove  a  small  number  of  special  cases  of  the  QAP  having  not 
only  polynomial-time  solutions,  but  solutions  that  can  be  expressed  in  closed-form  (so 
called  “easy”  cases).  See  Chapter  III  for  a  detailed  summary  of  their  findings. 
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2.  Sequencing  and  Scheduling 


The  body  of  literature  on  sequencing  and  scheduling  problems  is  rich, 
encompassing,  among  others,  the  areas  of  requirements  generation  (“open  shop” — where 
customer  demands  create  requirements  directly,  or  “closed  shop” — where  requirements 
are  in  the  form  of  inventory  replenishment  decisions),  processing  complexity  (one  stage 
or  multi  stage;  single  or  parallel  processors;  flow  shop  or  job  shop),  and  scheduling 
criteria  (schedule  cost,  schedule  performance).  Baker  [24]  provides  a  nice  introduction 
to  the  topic.  Graves  [25]  provides  a  fairly  comprehensive  review  of  the  production 
scheduling  literature  with  an  emphasis  on  these  three  areas.  Here  we  focus  on  the 
requirements  generation  portion  of  the  problem,  particularly  what  is  known  as  the 
sequential  assignment  problem. 

Derman,  et  al.  [26],  derive  a  recursively  defined  optimal  policy  for  assigning  n 
jobs  (arriving  in  sequential  order)  to  n  servers,  where  associated  with  each  job  j  is  a 
random  reward,  X . ,  and  associated  with  each  server  i  and  job  j  is  a  known  constant 

ptj ,  0<  ptj  <1  measuring  the  quality  of  the  server  ( ptj  =1  indicates  a  perfect  server). 
The  optimal  policy  assigns  the  n  jobs  to  the  n  servers  so  as  to  maximize  the  total  expected 

n 

reward.  Total  reward  is  defined  to  be  Albright  [27]  extends  the  results  in  [26] 

j= i 

to  account  for  discounted  rewards,  where  the  discount  is  a  function  of  time.  Agrawala,  et 
al.  [28],  derive  a  result  that  minimizes  expected  total  flow  time  (the  sum  of  job 
completion  times)  by  non-preemptively  sequencing  n  identical  jobs  (with  exponential 
service  times)  to  be  serviced  by  m  non-unifonn  processors.  Coffman,  et  al.  [29],  address 
the  problem  of  minimizing  expected  makespan  (maximum  job  completion  time)  by 
assigning  n  jobs  to  m  uniform  processors  with  exponential  service  times.  Optimal 
scheduling  rules  are  derived  for  two-  and  three-processor  systems.  Righter  [30]  derives 
results  that  minimize  a  weighted  function  of  job  completion  times  for  n  jobs  assigned  to 
m  processors  non-preemptively.  Service  time  depends  on  the  server  (not  the  job)  and  has 
increasing  failure  rate  (known  as  IFR).  Ross  [31]  states  an  optimal  policy  that  maximizes 
the  expected  total  reward  obtained  from  sequencing  n  jobs  for  assignment  to  a  single 
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processor.  Each  job  i  has  an  associated  independent  random  service  time,  X.,  and 
reward  a‘Ri  for  completion  at  time  t,  where  a  ( 0  <  a  <  I )  is  a  discount  factor. 
Voutsinas  and  Pappis  [32]  derive  a  similar  result  for  power  function  rewards  (i.e., 
functions  of  the  fonn  Vi  (l  )  =  K  lc,i ,  for  aj  <  0 ).  Some  of  the  analytical  results  we  derive 

in  Chapter  III  are  in  the  area  of  single-processor  sequential  assignment  problems,  and  are 
directly  related  to — and  in  some  cases  extensions  of — the  work  done  by  Righter  [30], 
Ross  [31],  and  Voutsinas  and  Pappis  [32]. 
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II.  MODEL  DEVELOPMENT 


In  this  section,  we  construct  two  models  that  relate  to  the  physical  system 
described  in  Chapter  I.  The  main  model  we  focus  on  is  a  queuing/job  sequencing  model 
(section  C);  a  companion  model — a  discrete  time  information  flow  model — is  introduced 
in  section  E.  We  begin  by  stating  the  necessary  definitions  and  assumptions. 

A.  DEFINITIONS 

Some  key  terms  are  defined  below. 

1.  Information  Package 

An  information  package  is  defined  as  a  discrete  unit  of  information  in  one  of 
many  forms  (e.g.,  a  single  text  message,  a  single  still  image,  or  a  single  video  clip). 
Information  packages  can  be  differentiated  by  processing  time  required,  size  (e.g., 
megabytes),  value  or  priority,  reliability  “decay”  rate,  source  (type  sensor,  configuration), 
and  other  relevant  factors.  Information  packages  are  represented  as  “jobs”  to  be  serviced 
(or  processed)  in  the  queuing/job  sequencing  model,  and  are  chosen  for  processing 
(individually)  based  on  their  value.  In  the  discrete  time  model,  information  packages  are 
represented  in  aggregate,  with  total  remaining  processing  time  (workload)  in  the  system 
the  metric  of  interest.  Some  examples  of  information  packages  include: 

•  A  two-minute  video  clip  received  from  a  high-altitude  unmanned  aerial 
vehicle  (UAV)  showing  suspected  militants  arriving  at  a  meeting  location 

•  A  30-second  voice  transmission  comprising  a  situation  report  on  the 
location  and  activity  of  a  threat  air  defense  artillery  battery 

•  A  three-megabyte  satellite  image  of  a  two-by-two  kilometer  zone  within  a 
city  of  interest 

Multiple  information  packages  can  describe  the  same  target.  Each  distinct  piece  of 
information,  whether  individual  reports  from  different  sources  or  multiple  reports  from 
the  same  source,  is  treated  as  an  individual  information  package. 
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2. 


Information  Package  Value 


We  assume  throughout  this  analysis  that  each  information  package  can  be 
assigned  a  value  function  that  describes  how  its  value  changes  over  time.  We  further 
assume  that  the  value  function  associated  with  an  information  package  is  a  non¬ 
increasing  function  of  time.  The  value  function  types  considered  throughout  this  analysis 
consist  of  the  following: 

•  Exponential:  Valj  (t)  = 

•  Linear:  Valj  (t)  =  //,  -  ayt 


Step:  Valj  (t)  = 


0  <t  <Tj 

otherwise 


Initial  value  and  the  rate  at  which  value  decays  for  each  information  package  are  defined 
by  the  function  parameters  j3,  a,  and  r .  For  example,  an  information  package 
containing  the  location  and  velocity  of  a  high  value,  time  sensitive  target  (e.g.,  a  sedan 
carrying  the  head  of  A1  Qaeda  in  Iraq)  could  have  a  value  function  that  is  exponential 
with  relatively  high  /?  (initial  value)  and  a  (value  decay  rate).  We  say  that  value  is 
realized  once  an  information  package  has  been  processed  (see  Figure  5  below). 


B.  ASSUMPTIONS 


We  make  the  following  assumptions  throughout  this  analysis.  Other  assumptions 
will  be  specified  as  appropriate. 

•  Realized  value  is  linearly  additive. 

•  Information  package  value  is  nonnegative. 

•  Information  package  arrival  and  processing  rates  are  constant  with  respect 
to  time  (i.e.,  there  is  no  “seasonality”  present  in  these  rates). 

•  Effects  on  the  system  due  to  bandwidth  constraints  are  not  considered. 

•  Service  will  not  be  preempted;  once  an  information  package  is  assigned  to 
a  processor,  it  will  not  exit  until  processing  is  complete. 

•  Processors  are  equally-capable  (i.e.,  “uniform”). 
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C.  QUEUING/JOB  SEQUENCING  MODEL  FOR  INFORMATION  VALUE 


1.  General  Model 

We  represent  the  battlefield  information  flow  system  as  a  multi-class  queuing 
system  with  general  inter-arrival  time  and  service  time  distributions  and  multiple 
servers/processors  (denoted  a  GIGIm  multi-class  queue),  with  the  following  defined 
characteristics: 

•  j  =  1, 2, . . . ,  J  information  package  (IP)  index. 

•  k=\,2,...,K  class1  index. 

•  Sj  service  time  length  for  IP  j  (random  variable). 

2 

•  X.  arrival  rate  for  IP  j. 

•  jtr  service  rate  for  IP  j? 

•  Aj  arrival  time  of  IP  j  (random  variable) 

•  Xj  the  time  IP  j  is  assigned  to  a  processor  (this  is  our  control).  We 
assume  that  processing  begins  immediately  upon  assignment. 

•  Dj  service  completion  (system  departure)  time  for  IP  j.  D/  is  function  of 
both  Xj  and  Sj  (specifically,  Dj  =  Xj  +  $ , ) 

•  Valj  (t)  value  of  IP  j  at  time  t 

Often  in  queuing  applications  we  assume  exponential  inter-arrival  and  service 
times;  in  this  case,  Xj  and  ///  represent  the  parameters  of  these  distributions.  The 

models  described  in  this  section,  however,  do  not  require  the  assumption  of  exponential 
inter-arrival  and  service  times;  we  therefore  consider  Xj  and  //  ,  as  rates  in  the  general 

sense. 


1  The  number  of  classes  in  the  physical  system  is  small  relative  to  the  number  of  information  packages 
present;  common  examples  of  information  package  classes  include  video,  still  imagery,  voice,  and  text. 

-  Arrival  and  service  distributions  (and,  hence,  rates)  are  thought  of  as  functions  of  IP  class  (thus  the 
system  has  at  most  K  unique  arrival  and  service  time  distributions);  however,  for  notational  simplicity,  we 
assign  arrival  and  service  rates  to  each  IP  (thus  the  subscript  j  rather  than  k). 
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Prior  to  entering  the  queue,  information  packages  (IPs)  are  assumed  to  undergo  a 
“pre-processing”  stage,  where  value  functions  and  class  are  assigned  to  each  information 
package;  see  Figure  4.  Value  functions  are  piecewise  continuous  functions  of  time,  and 
are  described  in  detail  in  Section  A. 


Figure  4  G/G/m  Multi-class  Queuing  System  Diagram 


Information  packages  arrive  into  the  queue  at  rate  Xj  and  wait  to  be  selected  by  a 

controller  for  processing  by  one  of  m  servers,  assumed  to  be  equally  capable.  We  assume 
that  infonnation  package  arrival  rates  far  outpace  the  system’s  processing  capacity — 
rendering  the  system  inherently  unstable  (queue  length  grows  without  bound,  along  with 
expected  wait  times) — which  mirrors  the  actual  command  and  control  system  under  study 
(which  we  will  refer  to  as  the  “physical  system”). 

The  role  of  the  controller  is  to  establish  a  dynamic  service  policy  that  results  in 
the  maximum  value  “realized”  by  the  system.  To  be  more  precise,  we  define  the  realized 
value  associated  with  package  j  as  follows: 
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V;  =  Vali(dj)  =  Valj(zj+sj) 


(2.1)3 


As  this  definition  implies,  realized  value  is  computed  upon  completion  of 
service — see  Figure  5.  The  departure  time  of  a  package  is  not  deterministic,  however,  so 
we  define  expected  realized  value  associated  with  package  j  in  (2.2). 


E [v; ]  =  E [Valj  (Dj )]  =  E[valj  (Zj  +  SJ ) 


(2.2) 


Figure  5  Measuring  Realized  Value 


Finally,  we  define  the  cumulative  expected  realized  value  obtained  from 
processing  a  set  of  Q  information  packages  as  follows: 


£[>']  =  £  tval,(D,)  + 

j= 1  7=1 


(2.3) 


3  The  lower-case  d-  and  S-  are  used  here  to  denote  realizations  of  the  random  variables  D.  and  S  ■ , 
respectively. 
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2.  Simplified  Model  and  Variations 


We  start  with  a  simplified  model  of  the  queuing  system  in  order  to  obtain 
analytical  results.  To  begin,  we  assume  a  GIGIm  queue  with  time  independent  arrival  and 
service  times  and  no  service  preemption.  We  also  assume  zero  time  delays  due  to  pre¬ 
processing  (assigning  information  package  class  and  value  function  prior  to  entrance  into 
the  queue),  server  assignment  decisions,  and  bandwidth  constraints.  Finally,  we  assume 
that  value  is  additive;  in  other  words,  value  realized  from  processing  various  information 
packages  can  be  linearly  aggregated. 

The  decision  process  starts  when  a  processor  becomes  available,  a  time  which  we 
will  call  td  (or  decision  time).  Let  Q  be  the  number  of  information  packages  present  in 

Oval  2  at  time  td  4  An  iterative  approach  is  employed  as  follows. 

•  Sequence  the  Q  information  packages  in  such  a  way  as  to  maximize 
cumulative  expected  realized  value  obtained,  assuming  all  Q  information 
packages  will  be  processed.  This  sequence  becomes  the  initial  processing 
“plan”  or  strategy,  and  may  be  updated  as  necessary. 

•  Begin  executing  this  processing  strategy  by  processing  the  first 
information  package  in  sequence. 

•  When  a  processor  next  becomes  available,  decide  to  either  continue  with 
the  original  strategy,  or  modify  this  strategy  by  updating  the  sequence. 
The  choice  to  update  the  sequence  is  made  if: 

0  New  information  packages  have  arrived,  or 

0  Actual  service  time  of  the  first  information  package  differs 
significantly  from  the  expected  service  time. 

•  If  updating  is  warranted,  redefine  Q  as  appropriate  and  re-sequence. 

•  This  new  sequence  becomes  the  “current”  plan  or  strategy,  and  execution 
is  begun  on  it  by  sending  the  first  package  in  sequence  to  the  available 
processor. 

•  Continue  iterating  in  this  fashion  until  stopping  criteria  are  met. 

In  this  approach,  we  update  (or  at  least  consider  updating)  the  current  strategy 
every  time  a  processor  becomes  available.  Of  course,  other  update  alternatives  exist; 

4  Q  may  instead  represent  a  subset  of  the  information  packages  present  in  Oval  2  at  time  t,  if  it  makes 
sense  operationally  to  omit  some  information  packages  from  consideration. 
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e.g.,  we  could  choose  to  abide  by  the  original  plan  for  a  set  amount  of  time  or  for  a  set 
number  of  information  packages  before  updating  the  plan/sequence.  We  choose  the 
updating  scheme  outlined  above  because  it  most  closely  matches  the  strategy  that  would 
be  employed  in  the  physical  system. 

We  now  state  the  optimal  control  problem  for  any  set  of  Q  information  packages. 
The  solution  to  this  problem  is  the  processing  sequence  that  maximizes  cumulative 
expected  realized  value. 

Let  Q  be  the  number  of  information  packages  in  the  set,  let  q  =  \,2,...,Q  be  the 
identifying  index  for  each  infonnation  package,5  and  let  g  =  1,  2,...,  m  represent  the 
servers  present  in  the  system.  Additionally,  define  the  following  terms: 

•  ng  number  of  IPs  assigned  to  server  g  (note  that  2>s=e) 

•  jc  .  e{l,2 IP  to  be  processed  hth  by  server  g  (this  is  the 
“decision  variable”) 

•  cpQ  the  set  of  all  possible  pennutations  of  the  sequence  (l,  2, . . . ,  Q ) 

Then,  the  optimal  control  problem  for  the  queuing/job  sequencing  model  is  as 
follows: 


max 

KaWe 


E[vr~ 


LL£ 

g=l  h= 1 


Val 


KJ 


(2.4) 


This  model  seeks  the  processing  sequence  of  information  packages  (represented 
byx  ,  )  that  results  in  the  maximum  cumulative  expected  realized  value,  and  will  be 

developed  in  detail  in  Chapter  III.  Figure  6  below  depicts  a  typical  information  package 
assignment  outcome. 


5  q  is  a  “temporary”  package  index;  each  package  still  maintains  its  permanent  index,  /. 
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As  we  shall  see  in  Chapter  III,  (2.4)  is,  in  general,  NP-hard.6  We  develop 
solutions  for  some  special  cases  in  Chapter  III,  and  resort  to  heuristic  strategies  for  more 
complex  cases  in  Chapter  IV. 

D.  DATA 

For  the  simulation  experiments  described  in  Chapter  IV,  we  require  reasonable 
estimates  of  the  distributional  parameters  for  the  random  variables  defined  in  the  previous 
section.  We  turn  to  a  set  of  experimental  data  obtained  during  a  recent  series  of  C2 
experiments  conducted  by  the  Defense  Advanced  Research  Projects  Agency  (DARPA)  to 
at  least  partially  satisfy  this  requirement.  See  [33]  for  details  on  this  experiment.  Table  1 
depicts  an  excerpt  of  the  data  collected  during  this  effort.  Each  row  represents  an 
information  package,  with  associated  information  about  the  target  (enemy  element), 
friendly  spotter,  arrival  times  to  Ovals  2  and  4  of  the  Dynamic  Model  of  Situated 
Cognition,  and  the  total  time  it  took  the  package  to  go  from  Oval  2  to  Oval  4. 


°  A  problem  IT  is  NP-complete  if  IT  e  NP  and,  for  all  other  problems  IT '  e  NP  ,  IT '  transforms 
to  IT  .  A  problem  IT  is  NP-hard  if  there  exists  some  NP-complete  problem  IT '  that  Turing-reduces  to  IT 
[16]. 
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Info 

Pkg# 

Enemy 

Element 

Friendly 

Source 

First 
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Time 

(Ai  i  ived 
Ov.il  2) 

Friendly 

Spotter 

Time  HTR 
Classified 

(Ai  i  ived 
Oval  4) 

Time  from 
02-04 
(min) 

12 

DRA1-2 

IntelDump 

13:27:08 

ARVCAU2 

15:27:39 

120.5 

13 

PUR5-5 

IntelDump 

13:27:08 

UAV2CAU21 

15:32:44 

125  6 

14 

GAR  3-1 4 

IntelDump 

13:27:08 

HHQ0U4 

15:41:30 

134.4 

15 

DAR6-2 

IntelDump 

13:27:08 

UAV43 

15:49:19 

142.2 

16 

DRAM  2 
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UAV42 

15:58:00 
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17 
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UAV42 
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18 
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19 
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UAV41 
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UAV41 
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6.0 

20 
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14:00:22 
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21 
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13:51:21 
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22 
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UAV2CAU22 

13:53:44 

0.2 

23 

DRA1-D5 

UAV2CAU22 

13:53:30 

UAV2CAU22 

13:54:17 

0.8 

24 

FT4-33 

UAV1 CAU1 2 

13:54:38 

UAV1 CAU1  2 

13:54:56 

0.3 

25 

SA1 3-6-1 

UAV2CAU22 

13:56:00 

UAV2CAU22 

13:56:17 

0.3 

26 

URA4-1 

UAV41 

13:57:05 

UAV41 

14:27:17 

30.2 

27 

URA8-2 

ARVCAU1 

13:57:29 

ARVCAU1 

13:57:36 

0.1 

Table  1  Excerpt  from  DARPA  C2  Experiment  Data 

(from  [34]) 


Data  in  the  column  titled  “First  Spotted  Time”  are  used  to  estimate  inter-arrival 
time  distributions,  and  data  in  the  column  titled  “Time  from  02-04”  are  used  to  estimate 
service  time  distributions.  See  [34]  for  the  assumptions  inherent  in  this  data  set  and  other 
details. 

Figure  7  below  depicts  a  histogram  of  inter-arrival  time  data  from  the  DARPA  C2 
experiment  along  with  the  best- fit  exponential  distribution  (with  parameter  X  =  0.63); 
Figure  8  shows  the  same  for  service  time  data  (with  parameter  ju  =  0.134  ).  We  conclude 
that  it  is  reasonable  to  model  inter-arrival  and  service  times  exponentially. 
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Density  «  Density 


Data 

7  Inter- arrival  Time  Data  with  Best- fit  Exponential  Distribution 


Figure  8  Service  Time  Data  with  Best-fit  Exponential  Distribution 
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E. 


DISCRETE  STOCHASTIC  INFORMATION  FLOW  MODEL 


We  now  develop  a  companion  model  to  the  queuing/job  sequencing  model  that 
represents  the  dynamics  of  information  volume  (or  workload)  in  the  system.  Figure  9 
below  illustrates  how  infonnation  passing  from  Ovals  2  to  3  of  the  Dynamic  Model  of 
Situated  Cognition  first  must  pass  through  an  information  processing  stage.  This  is  the 
specific  context  within  which  the  dynamical  system  model  of  information  flow  resides. 


Ovall 

9 


Sensor 

data 


Information 

Processing 


Processed 
data  j 

\ 

Oval 

3 

/ 

Figure  9 


Information  Processing  Layer  in  the  DMSC 


1.  Definitions  and  Assumptions 

To  develop  the  model,  we  first  require  some  definitions  and  assumptions: 

•  required  processing  time — a  pre-assigned  length  of  time  associated  with 
each  information  package  (e.g.,  a  particular  still  image  requires  10  time 
units  to  process). 

•  jc.  (k)  — the  number  of  information  packages  in  Oval  2  at  the  start  of  time 
period  k  for  which  the  remaining  processing  time  required  is  i  time  units. 

•  ui  (k  )  — the  number  of  information  packages  arriving  in  Oval  2  during 
time  period  k  which  require  a  totaV  of  i  time  units  to  process. 

•  At — time  step  size. 

•  V (k  ) — total  volume;  a  time-dependent  measure  of  the  workload  present 

in  the  system  (obtained  by  summing  the  remaining  processing  time  values 
of  each  information  package  present  in  the  system  (Oval  2)  for  each  time 
period  k). 


7  Note  the  difference  between  total  and  remaining  processing  time  required  in  the  definitions  of 
x^k)  and  ul  ( k ) .  When  an  information  package  enters  the  system,  it  arrives  with  an  assigned  total  required 
processing  time — and  the  clock  starts.  As  time  progresses  in  the  model,  the  required  processing  time  is 
decremented  appropriately — yielding  the  remaining  processing  time  required  tracked  by  xt(k) . 
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2. 


Model  Development 


We  now  define  the  Discrete  Information  Flow  Model  as  follows: 


x(k  + 1)  =  Ax(k")  +  Bu(k) 


(2.5) 


where  matrices  A  and  B  are  parameters  of  the  system,  with  B  initially  defined  as  the 
identity  matrix,  and  A  (known  as  the  “upshift”  matrix)  initially  defined  as  follows: 


A  = 


,  /  =  1 . . .  m,  j  =  1 . . .  n,  where  atJ  = 


jl,  ifj  =  /  +  l 

[0,  otherwise 


(2.6) 


The  vectors  x  and  u  contain  n  elements — corresponding  to  the  n  discrete 
processing  time  “bins”  of  size  At .  The  m  rows  of  A  correspond  to  the  m  time  steps  of  size 
At  represented  in  the  model  (i.e.,  in  •  A/  i s  the  total  time  modeled).8  A  is  a  matrix 
containing  all  zeros  except  for  the  off-diagonal  immediately  above  the  main  diagonal, 
which  contains  all  ones.  A  is  necessarily  square,  meaning  we  must  restrict  ourselves  to 
cases  where  m  =  n.  This  defines  the  completely  deterministic  (and  most  simplistic)  case, 
where  the  required  processing  time  for  each  information  package  in  the  system  is 
decremented  Ai  time  units  every  time  step,  and  the  number  of  packages  arriving  each 
time  step  (it)  is  fully  known. 

We  now  address  a  more  complex  case,  where  the  input  vector  u  is  a  function  of 
two  random  processes — one  each  for  information  package  arrival  and  processing  times. 
We  can  think  of  this  as  analogous  to  a  birth-death  process  [35],  where  information 
packages  arrive  (birth)  and  depart  (death)  from  Oval  2  at  known  (or  estimated)  rates. 


8  Processing  time  bin  size  and  time  step  size  may  be  chosen  so  as  to  not  be  equal;  in  general,  we  will 
keep  them  of  equal  size. 
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Define  the  following  random  variables: 

•  Yj  -  time  between  arrivals  of  packages  /-I  and/  (referred  to  as  inter¬ 
arrival  times). 

•  7 ’.  =  ]jT  7  -  arrival  time  of  package  j. 

9=1 

•  f)  -  number  of  packages  arriving  in  Oval  2  in  time  period  k. 

•  i? .  -  required  processing  time  for  package  j. 


Then, 


(\T~ 

j 

—  k 

At 

V 

7 

,  where 


x  =  a 
x^a 


(2.7) 


The  stochastic  version  of  the  discrete  infonnation  flow  model  now  becomes 


x(k  + 1)  =  Ax(k^  +  Bu(k} 


where 


At 


-/ 


j 


(2.8) 


Finally,  we  can  fonnally  define  total  volume  as 


V(k)  =  -£i-x,(n  (2.9) 

1=1 

3.  Model  Validation 

To  test  how  well  the  stochastic  version  of  the  discrete  infonnation  flow  model  fits 
the  actual  data,  a  simple  Monte  Carlo  simulation  is  developed  (see  Appendix  B  for 
MATLAB  code).  We  utilize  the  metric  V (k  )  (volume)  to  facilitate  this  comparison. 
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Each  run  of  the  simulation  proceeds  in  the  following  manner: 

•  Randomly  draw  inter-arrival  times  and  service  times  from  the  appropriate 
distributions9 

•  Choose  At  (At  =  2  min  for  this  assessment) 

•  Compute  u(k),  k  =  1  ...m  from  (2.7)  and  (2.8). 

•  Compute  x(k),  k  =  1  ...m  from  (2.5)  and  (2.6). 

•  Compute  V(k),  k  =  l...m  from  (2.9). 

Figure  10  displays  the  results  of  1000  runs  of  this  simulation  plotted  with  the 
original  experiment  data. 


Figure  10  Simulation  Results 


9  in  this  case,  we  used  Gamma  distributions  (with  parameters  estimated  from  the  DARPA  C2 
Experiment  data)  to  generate  inter-arrival  and  service  times.  It  turns  out  that  the  Gamma  distribution 
produced  a  slightly  better  fit  to  the  DARPA  data  than  did  the  exponential  distribution  for  both  inter-arrival 
and  service  times. 
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This  graph  plots  total  volume  (workload)  in  the  system  over  time  for  simulated 
results  and  actual  experiment  data.  The  solid  line  represents  the  mean  total  volume  for 
the  1000  simulation  runs;  dashed  lines  are  ±1.5  standard  deviations  from  the  mean.  The 
“x’s”  represent  the  original  DARPA  experiment  data.  From  this,  we  observe  that  the 
stochastic  model  consisting  of  (2.5)  and  (2.8)  is  a  reasonable  representation  of  the 
dynamics  and  random  processes  present  in  the  actual  data. 
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III.  ANALYSIS  AND  COMPLEXITY  RESULTS 


We  now  develop  the  model  (2.4)  in  some  detail.  All  cases  in  this  chapter  assume 
a  single  processor  ( m  =  1).  We  redefine  the  decision  variable  to  reflect  this  assumption 
as  follows:  x  is  the  information  package  to  be  processed  c/h.  All  results  presented  are  of 

the  “static”10  type;  i.e.,  we  assume  a  set  of  infonnation  packages  of  size  Q  and  seek  to 
process  them  in  the  sequence  that  results  in  maximum  realized  value  attained.  We 
assume  no  new  arrivals  and  no  service  preemption.  Relative  to  the  most  realistic  cases  of 
the  physical  system,  this  problem  is  very  much  oversimplified.  However,  this 
simplification  is  necessary  for  optimality  and  complexity  analysis,  which  provides 
fundamental  insights  for  the  overall  problem  and  supports  further  exploration  of  practical 
solutions. 


A.  STEP  VALUE  FUNCTIONS 

Consider  the  case  where  information  package  value  functions  take  on  the  form  of 
step  functions;  that  is,  consider  value  functions  of  the  form 


Valq  (0  =  pqH  (t  -  aq )  H  (aq  +  rq  - 1) 

\Pq  aq^t^aq  +  Tq 

0  otherwise 


(3.1) 


where  //  ( /  )  is  the  unit  (or  Heaviside)  step  function,  aq  is  the  arrival  time  of  information 
package  q  and  J3  ,  r  >  0 .  See  Figure  1 1 . 


10  See  Chapter  IV  for  details  on  static  and  dynamic  modeling. 


27 


Valq{i) 


a. 


a 


+  T 


x, 


t 


Figure  1 1  Step  Value  Function 


We  assume  that  a  <td  <a  +  r  ,  \/q.  Suppose  also  that  service  times  sq  are 
deterministic  and  that  m  =  1 .  Then  maximizing  V  is  equivalent  to 


(3.2) 


(3.3) 


(3.4) 


Expression  (3.3)  equates  to  maximizing  the  combined  realized  value  of  the 
information  packages  that  can  be  processed  prior  to  their  value  becoming  zero  (at  time 
t  =  aq  +  xq ),  while,  equivalently,  (3.4)  seeks  to  minimize  the  combined  realized  value  of 

information  packages  that  cannot  be  processed  “in  time.”  We  will  refer  to  (3.4)  as 
Problem  I. 
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Theorem  3.1:  Problem  /  is  NP-hard. 


Proof.  We  first  define  the  decision  problem  corresponding  to  the  optimization 
problem  /:  Let  K  >  0 .  Does  there  exist  a  permutation  (xx ,  x2 , . . . ,  xQ  ^  e  (pQ  such  that 


(3.5) 


We  will  refer  to  decision  problem  (3.5)  as  /’ . 

Consider  now  the  known  NP-complete  problem  “Sequencing  to  Minimize  Tardy 
Task  Weight”  (short  name  SS3)  in  [16],  defined  as  follows: 

Let  q  =  \,2,...,Q  tasks,  each  having  length  l(q)e  Z+ ,  weight  and  a 

deadline  d(q)eZ+.  Let  K  be  a  positive  integer.  Is  there  a  permutation 
(xr ,  x2 , . . . ,  xQ  ^  e  cpQ  such  that 


(3.6) 


Problem  SS3  (3.6)  is  clearly  a  special  case  of  problem  /' ,  where  l(q) ,  w(q) ,  and 
d(q)  (all  positive  integers)  are  analogous  to  sx  ,  J3q,  and  aq  +  tq ,  respectively.  The 

listed  quantities  for  Problem  /',  however,  are  not  restricted  to  integer  quantities.  It 
follows,  then,  that  problem  /’  is  NP-hard  (the  general  problem  is  at  least  as  hard  as  a 
special  case).  Furthermore,  we  conclude  that  Problem  /  is  at  least  as  hard11  as  Problem 
/  ’ ,  and,  hence,  NP-hard  as  well.  □ 


1 1  If  the  solution  to  Problem  /  is  known,  then  the  solution  to  Problem  /’  is  known  for  any  K.  The 
converse  is  not  true,  in  general. 
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Consider  now  the  case  where  information  package  value  functions  take  on  the 
form  of  “multi-step”  functions,  shown  in  Figure  12  below.  Corollary  3.2  follows 
immediately  from  the  fact  that  the  value  function  in  (3.1)  is  a  special  case  of  the  one  in 
Figure  12. 


Valq{t) 


aq  U  aq  +Tq 

t 

Figure  12  Multi-step  Value  Function 

Corollary  3.2:  The  optimal  control  problem  involving  information  packages  with 
multi-step  value  functions  is  NP-hard. 

B.  LINEAR  VALUE  FUNCTIONS 

We  now  present  a  set  of  results  specifying  the  optimal  processing  strategy  for  a 
single  server  ( m  =  1 ),  where  all  information  packages  have  piecewise  linear  value 
functions  of  the  fonn 


Valq{t)  =  H[t-aq)H[zq-t)\_Pq-aq{t-aq)\ 

0,  otherwise 


(3.7) 
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//  (/)  is  the  unit  (Heaviside)  step  function,  aq>  0  is  the  rate  at  which  value  decays,  and 
P  >  0  is  the  initial  value  for  each  information  package.  Note  that  Valq  (/)  =  0  whenever 

t<a  or  t  >  x  =  —  -  a  .  See  Figure  13  below. 

1  l  cc 

Ug 


Figure  13  Linear  Value  Function 

Based  on  Theorem  3.1  and  Corollary  3.2,  we  conjecture  that  the  optimal  control  problem 
involving  value  functions  of  the  type  specified  in  (3.7)  is  NP-hard.  However,  as  we  shall 
see  below,  this  problem  is  rendered  solvable  when  one  additional  condition  is  met. 

1.  Main  Results 

We  assume  that  the  Q  information  packages  have  all  arrived  by  time  td .  If  we 
also  assume  that  Q  is  chosen  such  that  all  Q  information  packages  can  be  processed  prior 
to  time  t  =  r  ,  \/q  e  {l,2,...,Q}  ,12  we  can  simplify  the  value  function  as  follows. 

Valq{t)  =  pq-aq[t-aq)  (3.8) 

Not  an  unreasonable  assumption,  it  turns  out.  This  assumption  can  hold  for  combinations  of  small 
Q,  small  a  ,  and  large  Pq.  We  will  see  cases  in  Chapter  IV  where  simulation  results  affirm  the 
plausibility  of  this  assumption. 
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Note  that,  without  loss  of  generality,  we  can  shift  the  coordinate  axes  such  that  td=  0 . 
With  fiq  redefined  appropriately,  the  value  function  for  each  infonnation  package  can  be 
simplified  as  follows. 


ValAt)  =  Pq~aqt 


(3-9) 


The  total  expected  realized  value  obtained  from  processing  all  Q  infonnation  packages  is 
then 


E\vr^  =  E 


q=\  V  h= 1  J 


which  is  the  quantity  we  wish  to  maximize.  It  follows  that 


q= 1  |_  h= 1 


Q 

.?=! 


tt«AA 

q= 1  /z=l 


(3.10) 


(3.11) 


(3.12) 


where  the  first  term  in  (3.12)  represents  the  total  initial  value  for  all  infonnation  packages 
and  is  constant  with  respect  to  x.  So,  maximizing  £’[Fr]  is  equivalent  to 


min 

{*i ,}^<Pq 


lix, £[X]- 


(3.13) 


This  is  the  optimal  control  problem.  We  now  present  two  results  based  on  (3.13): 
Lemma  3.3,  an  intermediate  result,  assumes  only  two  information  packages  to  be 
sequenced.  Lemma  3.4,  the  main  result  in  this  section,  assumes  any  number  of 
information  packages. 
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Lemma  3.3.  Let  Q  =  2  and  m  =  1 ,  and  suppose  that  Q  is  chosen  such  that  all  Q 
information  packages  can  be  processed  prior  to  time  t  =  r  ,  \/q  e  .  Let  the 

value  function  for  information  package  q  be  Valq  (t)  =  f3  -  aqt,  q  =  1, 2  .  Then  processing 
information  package  1  prior  to  2  is  optimal  if  and  only  if 


<2j  a2 


(3.14) 


Proof.  The  expected  realized  value  from  processing  infonnation  package  1  then 
2  is 


E 


£[A-«A  +  /?2-«2(s1  +  s2)] 


(3.15) 


=  /?, -alE[Sl]  +  p2 -a2E[Sx  +52] 

Processing  the  infonnation  packages  in  the  reverse  order  yields  the  following  expected 
realized  value. 


E 


E[p2-a2S2+pi-al(S,+S2)~\ 


=  P2  - a2E[S2 ]  +  /?!- ccxE [5)  +  S2 


(3.16) 


Then 


V, 


V). 


>E 


v, ; 


(24>. 


<=>  - axE [Sx \  +  P2~  oc2E[Sx  +S2]>  fi2- a2E [S2]  +  axE [5t  +  S2 


<=>  axE[S2 ]  >  a2E [5*, ] 


a,  a*, 

r  n  >  -  □ 
£[S,]  £[S2] 
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Lemma  3.4.  Let  Q>  2  and  suppose  the  remaining  assumptions  in  Lemma  3.3 
hold.  Then  processing  the  information  packages  in  decreasing  order  of  aqjE^Sq~^  is 
optimal. 


Proof.  Suppose  the  processing  sequence  ^qx,q2  ...,qj,qi,...,qQ^  is  optimal,  and 
that  aqp[Sq]>aq  E  Sq  .  Then  the  expected  realized  value  obtained  from 


processing  the  information  packages  in  this  sequence,  E 


V!' 


,  is  equal  to 


+£k-«,K+-+\+s«)' 


+  •••  +  £ 


+  •••  +  £ 


f><ra',(s«+-+s.,\ 
~a,g{s„  +,"+'st0) 


Similarly,  if  we  were  to  switch  the  processing  order  of  infonnation  packages  qj  and  qt 


(and  leave  all  others  unchanged)  the  expected  realized  value,  E 
be 


v; 


,?2,  .9, 


would 


+e[e-%{s„+-+s„+s,)_ 


+  ---  +  E 


+  ---  +  E 


p«-a,K+-+s«) 
~a,Q(s„  +  '"+'s®) 


The  optimality  assumption  yields, 
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E 

v;  J 

>E 

v;  . 

£[A,  -“«A,]+£  A,  _“«,(A  +  A)  +'"  +  £  ft.,.  iS,„  +"'  +  'S),  ) 


+  •••  +  £ 


•£[a-“,(a+-+s,as«)' 
£K-sa]+£[a-«*(a+5J' 
-E[p,ra<As„+-+s„+s,\ 


+  ---  +  E 


p1Q~a<,{s1'+'"+s<«). 

-•••+£[A-«,(A+"'+A) 

Av  ~“»c  I'5',:,  +  "'+ a) 


a  - (a+-+a)]+£ k  -  A  ( A  + " ' + A + A )' 


> 


A  (A  +'"+ A)]+£k  (A+-+A  +  A) 


-a  5 

■-a  S 

-E 

-a  S 

■-a  S 

-a  5 

A  <h 

a  aJ 

A  A 

A  A 

A  Aj 

> 


•~aS  1 

-£ 

5  — 

■-a  S 

L  Qi  <h 

<ii  <ii  j 

qi  q\ 

A  A 

Vi  qi 

V» 


which  is  a  contradiction.  □ 

The  proof  above  uses  a  technique  known  as  “adjacent  pairwise  interchange 
methods,”  found  in  [24]  and  [31],  among  others. 

2.  Quadratic  Assignment  Problem  Formulation  and  Result 

We  can  reformulate  the  optimal  control  problem  for  linear  value  functions,  (3.13), 
as  an  equivalent  quadratic  assignment  problem  (QAP)  as  follows.  Let  y..  be  an  indicator 

variable  that  equals  one  if  information  package  i  is  processed  /h  and  zero  otherwise. 
From  (3.13)  we  have 
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q—\  h= 1  h= 1  q=h 

=  E[s*](a*l+a*I  +"-  +  %) 
+  E[SXi\aXi+aXi+---  +  aXg) 

+  ••• 

+  E  Sx  (a  +a  ) 

e-i  J  V  xq- i  xe  / 

+£Mk) 
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CC ,  +  0T,  +  O',  H - +  Of 


£[S,]  £[S2]  £[S,]  -  £[sj 


-*1  *2  -*3 


or  +  0;  H - 1- or 

x2  x3  x 


O',  H - bor 


or. 


where  Y  =  v„ 

lJ  - 


E[S2]  •••  or2  ...  aQ]  YL 


where  L  = 


£[S,]  £[S2]  -  £[sj 


YLYr 


a , 


or,, 


or. 


nr 


1  (jll^^l]"1"  J2,^2]  "+  Jgl^j^g])  (  Jl  I  “*■  ^  Jig)  H  I"  OTg  (  Jgi  +  ■ 


"(  Jl2^  [*^1  ]  y 22^  [^2  ]  “*■  ^  Jg2'^'[‘^g])  ^1  (  Jl2  "*■  ^  Jig  )  ^  ^  a0  (  Jg2  ' 


'  (  Jlg^  [Si  ]  +  J2g^’[^2  ]  +  •••  +  Jgg^  [‘S’g  ])  a\  (  Jig  )  H  h  °Q  (  Jgg  ) 


Let  G  =  \J\T  =  [g(i  ]  = 


£j* 


l=k 


.  Then  (3 . 1 7)  is  equivalent  to 


(3.17) 


0  •••  0 

1  ■■.  : 

;  ■•.  o 

1  •••  1 


•+Jgg)_ 

••+Jgg) 


37 


I 


k= 1 


f  Q  \f  Q 

a,glt 

\i= 1  J\‘= 1 


X 

7 


£=1  i=l  y=l 

see  q 

=^LzLaJE[si]y*2Lyji 

k-\  i=l  y=l  /=& 

0  G  g  G 

=YY,YYaaAs]y  iky.ii 

i=l  7=1  k=l  l=k 


Q  Q  Q  Q 

=  ZIEEc9toi/’  where  c,y=«/[^]  and  dkl 

i=l  y=l  i=l  /=1 

Hence,  an  equivalent  optimal  control  problem  (to  (3.13))  is 

e  e  e  e 

minimize  ££££  cijdk,y,kyji 

i=l  7=1  k= 1  /= 1 
G 

subject  to  =1,  V/ 


E^=1» v/ 

;=i 

yv  e{o,i},  v/,7 


r  i,  if&</ 

[O,  otherwise 


(3.18) 


(3.19) 


This  is  a  special  case  of  the  well-known  quadratic  assignment  problem  (QAP) 
[14].  The  general  QAP  is  known  to  be  NP-complete  [15,  16].  We  show  here  that  not 
only  is  (3.19)  not  NP-complete,  its  solution  can  be  expressed  in  closed- form. 

Theorem  3.5.  Consider  the  quadratic  assignment  problem  (3.19),  where 
C  =  [c,]  =  \jXjE  [.S'-  ]J  .  If  the  Q  information  packages  are  indexed  such  that 

«1  >  «2  >  >  aQ 

£[S,]  £[S,]  £[S„]’ 
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(which  results  in  no  loss  of  generality)  then  the  optimal  solution  is  Y*  =  [ ”y*.  J  =  I . 

The  proof  follows  immediately  from  the  equivalence  of  QAP  formulation  (3.19) 
and  that  of  (3.13).  □ 

Theorem  3.5  is  significant  in  that  it  proves  the  existence  of  an  analytical  solution 
for  a  case  of  the  quadratic  assignment  problem  not  previously  known  to  have  an 
analytical  solution.  A  small  number  of  special  cases  have  previously  been  shown  to  have 
analytical  solutions  (see  [18-21,  23]),  but  (3.19)  is  not  among  them.  In  general,  these 
special  cases  are  defined  in  tenns  of  the  matrices  C  =  [cv]  and  D  =  \dkl],  and  some  well- 
known  ones  are  summarized  as  follows. 


1 .  If  either  of  the  following  conditions  holds,  then  the  identity  pennutation  is 
optimal  [20]: 


c,s  ^  cjs >  csi  ^  csj »  d„  >  djs ,  and  d„.  >  dsi,  V\<i<j<Q  and  \<s<Q  (3 .20) 


cu  —  •  •  •  —  Cqq,  dj ,  >  . . .  >  dQQ ,  cfa  +  csk  c,,  +  csi,  and  dst  dts 
V  1  <i  <k<Q  and  1  <s,t<Q 


(3.21) 


2.  If  C  is  monotone  anti-Monge  (i.e.,  -C  is  Monge)13  and  D  is  a  symmetric  Toeplitz 
matrix  generated  by  a  benevolent  function,  then  the  permutation 

n  =  (1,3,5, 7, 9,. ..,10,8,6,4,2)  is  optimal  [17,  19,  21].  An  nxn  matrix 

D  =  [</..]  is  a  Toeplitz  matrix  if  there  exists  a  function 
f  :{-n  +  \,...,n-l)  — »Msuch  that  dtj  =  f(i-j)  for  1  <  i,  j  <  Q .  The  Toeplitz 
matrix  D  is  said  to  be  generated  by  the  function  f.  A  function 
f  :[-n  +  \,...,n- 1}  — »  M  is  called  benevolent  if  it  fulfills  the  following  three 


properties:  1)  /(-/)  =  /(/),  V  l<i<Q-l;  2)  /(/)</(/  +  !),  VI  <i< 


Q 

2 


-1. 


and  3)  f(i)  <  V  1  <  i  < 


Q 

2 


-1  [19]- 


3.  If  C  is  a  monotone  anti-Monge  matrix  and  the  elements  of  D  satisfy  the 
inequalities  below,  then  n  is  the  optimal  permutation  [23]: 


I3  See  (3.36)  for  the  definition  of  a  Monge  matrix. 

39 


(3.22) 


^Q+l-i,Q+\-j  di,Q+\-j  ~  d Q+\-i,Q+\-j  d Q+i-ij  —  0,  V  1  <  i  ^  j  < 


Q 

2 


dQ+l-i  Q+X-j  +  dQ+l_tj  di  Q+l_j  dtj  <  0,  diQ+]_j  dQ+1_Uj  <0,  VI  <i,j< 


Q 

2 


(3.23) 


<e+2-y  -dj  >  0,  d^yj-dj,  >  0,  V  1  <i< 


Q+ 1 


,  2 <i< 


Q  + 1 


i^j  (3.24) 


d o+2-i,Q+2-j  +  dQ+2-ij  dy  di  Q+2_j  >  0,  di  Q+2_j  dQ+2_tj  >0,  V  2  <i,j< 


Q 

2 


(3.25) 


It  can  be  easily  shown  that  matrices  C  and  D  defined  in  (3.18)  and  (3.19)  do  not 
satisfy  any  of  cases  1-3  above  (namely  because  C  is  not  anti-Monge).  If  we  assume  that 
the  service  times  (S',.)  are  identically  distributed,  we  could  construct  C  so  as  to  be  anti- 
Monge;  but  D  would  still  not  meet  any  of  the  conditions  outlined  above. 

3.  Linear  Assignment  Problem  Formulation  and  Result 

The  result  outlined  below  represents  a  special  case  of  Theorem  3.5  which, 
ordinarily,  would  not  warrant  special  consideration.  However,  as  an  application  example 
of  the  rich  literature  on  Monge  matrices14  and  the  linear  assignment  problem  (see  [14]), 
we  deem  it  worthy  of  inclusion  here.  This  result  specifies  the  optimal  processing  strategy 
for  a  single  server,  with  all  information  packages  having  linear  value  functions  of  the 
form 


Valq{t)  =  H(t-aq)H(rq-t)\pq-a(t-aq) 


pq-a[t-dq),  aq<t<zq 


0, 


otherwise 


(3.26) 


14  An  m  x  n  matrix  C  is  called  Monge  if  it  satisfies  the  so-called  Monge  property 
Cy  +  crs  <  cis  +  crj  for  all  1  <i<r  <  m,  1  <  j  <s  <n  [37,  38], 
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where  H(t )  is  the  unit  step  function,  a  >()  is  the  rate  at  which  value  decays  (assumed 
equal  for  each  information  package  in  this  case),  and  ft  >  0  is  the  (unique)  initial  value 
for  each  information  package. 

Corollary  3.6:  Suppose  at  time  td ,  Q  information  packages  are  chosen  for 

processing,  all  having  linear  value  functions  of  the  form  given  in  (3.26).  Suppose  also 
that  m  =  1  and  that  Q  is  chosen  such  that  all  Q  information  packages  can  be  processed 
prior  to  time  t  =  r  ,  \/q  e  {l,2,...,(9}  .  Then  processing  the  Q  information  packages  in 
order  of  increasing  expected  service  time  is  the  optimal  policy. 

Proof:  By  assumption,  the  Q  information  packages  have  all  arrived  by  time  td 
and  will  all  be  processed  by  time  r  ,  so  we  can  simplify  the  value  function  as  follows: 

Valq(t)  =  Pq-a(t-aq)  (3.27) 

Note  that  we  can  again  shift  the  coordinate  axes  such  that  td  =  0 .  With  f3q 

redefined  appropriately,  the  value  function  for  each  information  package  can  be 
simplified  as  follows. 


Valq(t )  =  J3q-  at  (3.28) 

Define  xh  e  {l,  2, . . . ,  Qj  to  be  the  information  package  processed  hth  (a 

modification  of  our  original  decision  variable  to  account  for  the  presence  of  only  a  single 
server).  Then  the  total  expected  realized  value  obtained  from  processing  all  Q 
information  packages  is 


<7=1  |_  h=\ 


which  is  the  quantity  we  wish  to  maximize.  From  (3.29)  it  follows  that 


(3.29) 
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(3.30) 


£M=Ia,-“I££IX 

9=1 


Q  q 

II 

g=l  /z=l 


=  IA-“II£K]  (3.30 


where  the  first  term  in  (3.3 1)  represents  the  total  initial  value  for  all  infonnation  packages 
(value  at  time  t  =  td)  and  is  constant  with  respect  to  x.  So,  maximizing  £'[l/r J  is 
equivalent  to 


K1 


e  q=\  h= l 


(3.32) 


We  can  transform  (3.32)  into  a  more  recognizable  fonn  by  introducing  a  slightly 
different  decision  variable.  For  i,j  =  \,2,...,Q,  let  y..  be  an  indicator  variable  that 

equals  one  if  information  package  i  is  processed /h  and  zero  otherwise.  We  have 


£l>£[XHe  q- i  •••  i]y 


where  Y  =  [ ~ ytj  j .  Therefore,  the  optimal  control  problem  can  be  written: 

Q  Q 

minimize  IIw 

f=l  7=1 
O 

subject  to  Sky=1’V7  (3-33) 

i=l 

Zk,=l,  Vi 

7=1 

y..  e{(U},  Vi,  j 


a£[Sj] 

a£[S,] 

a£[st)] 


42 


where 


and 


Cy  =a(Q-i  +  l)E[Sj 


c=h]  = 


a 


QE[S ']  QE[S2 ] 

(e-ijsft]  (q-i)e[s2] 


QE[s0 ] 

(e-i)£[s0] 


£[S2] 


(3.34) 


(3.35) 


This  is  the  classical  linear  assignment  problem  (LAP),  with  cost  matrix  C  defined 
in  (3.35).  While  polynomial-time  algorithms  exist  for  solving  the  general  LAP  (see  [36] 
for  a  recent  survey),  there  are  certain  special  cases  whereby  an  analytical  solution  may  be 
obtained  (see  [18,  37,  38]). 

Assuming  service  time  distributions  are  known,  we  can  index  the  Q  information 
packages  (without  loss  of  generality)  such  that  E [5]]  <  E [S,]  <  •••  <  E ,  thereby 

causing  the  matrix  C  to  possess  special  properties.  The  tenns  in  each  row  are 
nondecreasing  from  left  to  right  due  to  our  indexing  of  the  information  packages,  and  the 
tenns  in  each  column  are  strictly  decreasing  from  top  to  bottom. 

Hoffman  [38]  and  Burkard  et  al.  [37]  define  the  following:  an  m  x  n  matrix  C  is 
called  Monge  if  it  satisfies  the  so-called  Monge  property 

Cy  +  crs  <  cis  +  c  .  for  all  1  <  i  <  r  <  m,  1  <  j  <s  <n.  (3.36) 

Lemma  3.7  ([18,  37]):  a  matrix  C  is  Monge  if  and  only  if  its  elements  satisfy  the 
following  inequality: 

cij  +  L+ij+i  +cmj  fora11  (3.37) 
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Lemma  3.8  [18]:  if  the  cost  matrix  of  a  linear  assignment  problem  satisfies  the 


Monge  property,  then  the  optimal  solution  is  the  “identity  permutation”  (i.e., 

v=[v, ;]=n 


We  now  show  that  the  matrix  C  in  (3.35)  is  Monge. 


Cij  Ci+l,j+\  Ci,j+ 1  Ci+i,j 

=(e-i+i)ir[s/]+(e-(i+i)+i)£[s,,1]-(e-i+i)£[sM]-(e-(<+i)+i)£['sJ] 

=£[s,](e-i+i-(e-(i-+i)+i))+£[y.,](e-(<+i)+i-(e-<+i)) 

=  £[^]-£[S,„] 

<0 

Hence,  C  is  Monge  by  Lemma  3.7.  By  Lemma  3.8,  the  optimal  solution  to  (3.33) 
is  Y*  =  | ”y*.  ]  =  I .  The  interpretation  of  this  solution  is  to  process  the  Q  information 
packages  in  order  of  increasing  expected  service  time,  and  this  completes  the  proof.  □ 

C.  EXPONENTIAL  VALUE  FUNCTIONS 

Suppose  all  information  packages  have  value  functions  of  the  form 

ralt(i)  =  pte*  (3.38) 

Suppose  that  Q  =  2  and  m  =  1.  Suppose  also  that  Sq  ~  Exponential  ( juq  j .  Then,  for 

(hj)e<P2, 
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V, 


u. 


=  E[j8ie~a‘s‘]  + E 


fifi 


aAs*+s.  ) 


Then, 


oo  oo  co 

=  fii/d.^e~a,Sie~mdsi  +  Pj/LijUi J  J e  a,'s'+Sl’e~me  MjSj dsids] 


o  o 


■+ 


PiW, 


a,+E,  (aj+fi^ccj+tij) 


V, 


(1,2). 


>E 


V, 


(2.1). 


.  ^  P\E\  _j_  PiEiEi 


>  A/t  ,  A/lm 


«i+M  (a2+//1)(«2+//2)  «2+/L  («!+//, )(«!  +  m) 


<=> 


Aa 


(2i  +  //j 


Mi 


1  — 

V  (2j  +  //2 


> 


Pith 


f 


a2  +  El 


Ei 


1- 

V  ^2  +  E\  j 


P\E\a\ 


> 


PiEE-^  2 


(«!  +//1)(<21  +//2)  (a2  +//2)(«2  +//j) 


(3.39) 


Unfortunately,  this  expression  does  not  lend  itself  to  a  distinct  ordering  as  was  the 
case  for  linear  value  functions. 

Let  us  now  consider  the  case  where  Q>  2.  It  is  easily  shown  that  for 
(i,j,k)  e  <p3. 


VE 


j*) 


PiEi 


-  + 


PjEtEi 


-  + 


PkEiEjEk 


J  «, 


+  Ei  (aj  +  Ei  )(oCj+  Ej )  (a* ■  +  Ei  )(«*•+  Ej )  K  +  Ek ) 


For  arbitrary  0  and  for  (c[l,q2,...,qQj  e  (pQ , 


E 


v: 


Wl-?2 


. ?e). 


AAi  |  I  I  •  A _ 

a,  K a,)K  a)  (%  a)Kb  a  )•••(% a2) 
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(3.40) 


Q  m 


=z/',n 

m= 1  n= 1 


+V% 


From  (3.39)  and  (3.40)  we  present  the  following  results  for  Q  >  2  .  Lemma  3.9  applies  a 
known  result  to  the  problem  at  hand;  to  our  knowledge,  the  results  presented  in  Lemmas 
3.10  and  3.11  are  new. 


Lemma  3.9:  If  aq  =a,  \/q  then  the  optimal  policy  is  to  process  information 
packages  in  decreasing  order  of  PqHq  (a  version  of  this  result  is  presented  in  [31]). 

Lemma  3.10:  If  //  =  //,  Mq  (i.e.,  service  times  are  identically  distributed)  then 


the  optimal  policy  is  to  process  information  packages  in  decreasing  order  of 


Pqaq 


K+//)2 


Suppose  now  that  Sq  ~  Uniform(a?,&?),  for  0  <aq  <bq.  Then,  for  (z,y)  e  <p2. 


v; 


j). 


E[^e-a‘Si]  +  E 


’Aw. 


=4 


r  i  a 

V  bi  -  ai  J 


b:  b: 


dSi+Pj\\' 


(  i  A 


V  b,  ~  at  j 


f  1  1 
-  ds  ds , 


_  a.a,  cCjb: 

e  1 1  -e  11 


)  /I  (e  "  "  -  e~ajbj )  [ea‘at  -  e  ah ) 


a. 


ibi  ai)  {aJ'f(bt-ai){bJ-aj) 


After  some  algebra,  one  can  show  that  E 


V, 


<uy 


>E 


v, ; 


(2,1) , 


if  and  only  if 


/?,  (e-a'a'  -e~aA ) 
ax  (hj  -  a, ) 


g~ala2  _ g~a\b2 

a\  ip2  ~a2) 


> 


a2(b2-a2) 


£~a2a\  _ g~a2b\ 

a2  (hj  -a,) 


(3.41) 


Lemma  3.11  follows  directly  from  (3.41).  The  proof,  not  shown,  is  constructed  exactly 
as  that  for  Lemma  3.4,  employing  an  adjacent  pairwise  interchange  argument. 
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Lemma  3.11:  Suppose  that  S  ~  Unifomi(a,h)  (i.e.,  service  times  are  identically 

distributed)  and  each  information  package  has  a  value  function  of  the  form  (3.38).  Then 
the  optimal  policy  is  to  process  information  packages  in  decreasing  order  of 


D.  AN  UPPER  BOUND  FOR  MIXED  VALUE  FUNCTION  TYPES 

Let  Val,  (/  )  be  the  (arbitrary)  value  function  for  information  package 
i,  i  =  l,2,...Q ,  let  in  =  1  server,  and  define  the  permutation  mapping  (p{  j)  to  be  the 
information  package  processed  /\  Then,  from  (2.4),  the  optimal  control  problem  is 


maximize 

<p 


E\_Vr~ 


O 

I* 

j= 1 


Val 


<pU) 


(A/))' 


(3.42) 


which  is  equivalent  to 


maximize 

<p 


e[v~ 


O 

’I* 

j= 1 


Val 


<pU) 


+  S 


<p(J ) ) 


(3.43) 


This  problem  is  conjectured  to  be  NP-hard. 

However,  if  we  let  si  be  the  (deterministic)  service  time  for  information  package 

i,  then  we  can  formulate  an  upper  bound  for  the  solution  to  (3.43) — which  is  to  say,  an 
upper  bound  on  maximum  cumulative  realized  value.  With  detenninistic  service  times, 
the  optimal  control  problem  becomes 


maximize 

v 


(3.44) 
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Let  (p*{j)  denote  the  permutation  (or  sequence)  that  solves  (3.44),  which  is  the 
sequence  that  produces  the  optimal  solution  to  (3.44),  denoted  V*.  Define 
Val  (t) 

T.  (t)  = - — ,  a  scaled  version  of  the  value  function  for  information  package  i,  and  let 

si 

=  V(i)  +  V(2)  "I —  +  s<p*(i)  /Ul  scrv'cc  completion  time).  Then  an  equivalent 

expression  for  V  *  is 

f*=IvmW'T  <3-45> 

j= 1 

The  geometric  interpretation  of  (3.45)  is  the  sum  of  rectangles  with  width 

and  height  Tpt( /( (^ )  •  Figure  14  below  shows  an  example  of  this  for  Q  =  4.  The  curves 

in  the  figure  (exponential  in  this  case,  but  they  could  be  any  value  function  type) 
represent  the  functions  Tf  (t) ,  i  =  1, 2, . . . ,  4 . 


I'Vwb  S(p{2)  *b  ‘Vw 

Figure  14  Computing  V*  Geometrically 
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Define  the  envelope  function,  Tenv  (7)  =  max  (7) ,T2  (7), . . T0  (7)j .  This  brings 
us  to  the  following  proposition  regarding  an  upper  bound  on  Ff 

'e 

Proposition  3.12: 

0 

Figure  15  below  gives  the  geometric  argument  for  this  proposition,  where  the 
envelope  function  Tenv  (7)  is  shown  in  black. 


IV)  b — V)  +  V)  *b  V) 


Figure  15  Comparison  of  V*  to  J  Tenv  (t^dt 


Here  we  have  constructed  the  function  Tenv  (  t)  from  the  functions  T.  (7)  = 


Valt  (7) 


and  show  how  it  can  be  used  to  construct  an  upper  bound  on  V  * .  We  will  see  in  the  next 
chapter  (Section  E)  that  the  function  7)  (7)  itself  yields  a  good  lower  bound  on  Ff 
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IV.  INVESTIGATION  OF  CONTROL  STRATEGIES 


A.  INTRODUCTION 

In  Chapter  III,  we  introduced  some  provably  optimal  controlled  queuing  strategies 
for  certain  special  cases  of  the  system  in  question.  We  also  noted  that,  even  for  several 
seemingly  oversimplified  cases,  the  optimal  control  problem  is  either  provably  NP-hard, 
or  conjectured  to  be  so.  To  address  more  complex  cases  (i.e.,  cases  that  more  closely 
mirror  those  present  in  the  physical  system),  we  depart  from  the  search  for  optimal 
control  policies  and  introduce  a  set  of  heuristic  control  strategies.  We  design  two 
simulation  experiments  to  address  two  related,  yet  distinct,  problems:  static  scheduling 
of  a  set  of  infonnation  packages,  and  dynamic  scheduling  with  information  packages 
arriving  in  real  time.  The  static  scheduling  problem  is  actually  a  sub-problem  of  the 
dynamic  scheduling  problem;  see  Figure  16. 

In  the  first  experiment,  we  assume  that  a  finite  set  of  information  packages  are 
available  at  time  zero  for  scheduling,  and  that  no  new  information  packages  arrive.  This 
experiment,  known  hereinafter  as  the  static  scheduling  simulation,  is  used  to  measure  the 
affect  of  sequencing  these  infonnation  packages  according  to  various  strategies.  The 
purpose  of  this  static  simulation  is  to  narrow  down  the  set  of  heuristic  strategy  candidates 
to  be  considered  in  the  second  simulation  experiment  by  measuring  their  performance  (in 
terms  of  the  metric  defined  in  (4.1))  and  robustness  to  variations  in  the  experiment 
factors.  This  second  experiment,  known  as  the  dynamic  scheduling  simulation,  is  meant 
to  be  a  credible  emulation  of  the  physical  system:  random  arrival  and  service  times;  an 
“unstable”  queue  (i.e.,  the  number  of  information  packages  awaiting  service  grows  faster 
than  they  can  be  processed);  and  a  realistic  mix  of  value  function  types.  Its  purpose  is  to 
compare  promising  sequencing  strategies  to  each  other  and,  ultimately,  to  a 
representation  of  the  “status  quo” — how  these  processing  decisions  are  currently  made. 

Figure  16  below,  depicts  the  basic  process  governing  the  physical  system  and  the 
dynamic  scheduling  simulation.  Upon  the  availability  of  a  processor  (at  time  td  ),  a  set  of 
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size  Q  of  the  information  packages  in  Oval  2  is  selected  for  consideration  for  assignment 
to  the  available  processor.  For  our  purposes,  this  set  is  considered  to  be  the  entire 
contents  of  Oval  2.  The  Q  information  packages  are  then  sequenced  either  optimally  (for 
small  Q )  or  according  to  some  alternative  sequencing  strategy.  This  sequencing  step 
defines  the  static  scheduling  problem,  which  can  be  thought  of  as  a  sub-problem  of  the 
dynamic  scheduling  problem.  The  first  information  package  in  the  sequence  is  assigned 
to  the  available  processor,  and  the  process  repeats  when  the  next  processor  becomes 
available.  Newly  arriving  information  packages  are  considered  for  subsequent  processor 
assignment.  The  process  repeats  itself  until  some  pre-defined  stopping  criteria  are  met. 


Select  Q  IPs 
from  Oval  2 
(new  arrivals  are 
considered) 


Sequence  IPs 
according  to 
desired 
strategy 


Figure  16  Physical  System  Process  Represented  in  the  Dynamic  Simulation 


B.  HEURISTIC  CANDIDATES 

The  construct  of  each  simulation  described  in  this  chapter  centers  on  the  decision 
to  assign  a  specific  information  package  to  the  next  available  processor,  and  the  strategy 
to  employ  when  making  this  decision.  Suppose  a  processor  becomes  available  at  time 
td,  d  =  1,2,..., •  The  following  heuristics  define  some  potential  assignment  strategies; 
in  each  case,  a  unique  processing  sequence  of  the  Q  information  packages  is  achieved. 
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HI:  Assign  infonnation  packages  to  processor  in  descending  order  of 

£[>J 

H2:  Assign  information  packages  to  processor  in  descending  order  of  ft 
(initial  value). 

H3:  For  each  assignment  decision  d,  choose  the  information  package  with 
the  highest  current  value;  i.e.,  the  highest  value  of  Valq  (td ) . 


H4:  For  each  assignment  decision  d,  choose  the  information  package  with 
the  highest  value  at  the  end  of  expected  service  time;  i.e.,  the  highest  value 


)■ 


H5:  For  each  assignment  decision  d,  choose  the  information  package  with 


the  highest  value  for  the  expression - 


Valq[E\pq\xq  =td 


H6:  Assign  infonnation  packages  to  processor  in  descending  order  of 
linear  combinations  of  terms  as  follows. 


0  H6E:  r/ 
only). 

0  H6L:  a 


aA  ,1-7 


E[Sq]  P 


0  <  ij  <  1  (for  exponential  value  functions 


V  ,  1-7 


CA+  A 


0  <  rj  <  1  (for  linear  value  functions  only). 


»  J 


LIFO:  Last  in,  first  out. 

FIFO:  First  in,  first  out. 

Random:  randomly  (uniformly)  select  the  next  information  package  to 
send  to  processor. 


Some  notes  and  further  explanation:  1)  HI  is  obtained  directly  from  the  optimal 
solution  for  the  linear  value  function  case  (see  Chapter  III,  Section  B);  2)  LIFO  and  FIFO 
are  only  suitable  for  study  in  the  dynamic  simulation  case,  since  the  static  simulation 
does  not  consider  the  arrival  time  of  each  information  package;  3)  H6E  and  H6L  are  not 
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suitable  for  study  in  simulations  where  value  functions  are  of  mixed  type;  4)  H2-H4 
represent  varying  types  of  “greedy”  strategies,  ranging  from  the  “naive”  H2  to  the  more 
sophisticated  H4. 

C.  STATIC  SCHEDULING  SIMULATION 

The  static  scheduling  simulation  incorporates  the  following  assumptions:  single 
processor;  Q  information  packages  that  have  all  arrived  by  the  start  of  the  simulation;  no 
new  arrivals;  various  mixes  of  value  function  types;  general  processing  time 
distributions;  and  no  service  preemption.  We  define  the  metric  V'  ,  a  surrogate  for 
expected  total  realized  value,  as  follows: 


(4.1) 


This  metric  is  a  function  of  r .  (the  time  information  package  j  is  assigned  to  a 
processor),  which  is  the  control.  It  is  preferable  to  a  metric  that  incorporates  a  rigorous 

Q 

notion  of  expected  value — say,  ValJ  \  xj  +  S  )  — because  it  can  be  applied  directly 

j= i 

regardless  of  value  function  type  and  service  time  distribution.  The  simulation  is  static  in 
the  sense  that  it  represents  a  “snapshot”  of  the  physical  system  (or  the  dynamic 
scheduling  problem)  at  a  moment  in  time — in  particular,  the  time  a  processor  becomes 
available — and  replicates  various  decision  strategies  regarding  which  information 
package  to  send  to  the  available  processor.  See  Appendix  A  for  the  static  simulation 
MATLAB  code. 


The  static  scheduling  simulation  is  not  meant  to  replicate  the  physical  system. 
The  purpose  of  the  static  simulation  is  to  provide  a  mechanism  with  which  to  evaluate 
service  policy  heuristic  candidates  by  computing  V  for  the  information  package 
processing  sequence  suggested  by  each  heuristic.  The  result  will  be  a  set  of  promising 
heuristics  to  use  in  the  dynamic  simulation. 
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1.  Comparison  Study 


The  first  component  of  the  static  simulation  involves  a  comparison  of  V'  for  the 
processing  sequence  suggested  by  each  heuristic  with  that  for  an  optimal  sequence.  Due 
to  the  combinatorial  nature  of  the  optimization  problem — i.e.,  determining  which  of  the 
Q\  possible  sequences  maximizes  V' — optimal  solutions  can  only  be  found  for 
relatively  small  Q.  We  begin  by  examining  randomly  generated  sets  of  information 
packages  that  vary  according  to  the  following  parameters: 

•  Value  function  type:  exponential  or  linear 

•  Value  function  parameters:  a  ~  Uniform(0.1, 1)  and  /?  ~  Uniform(5,10) , 
where  Valq  (/)  =  f:Sq  -aqt  and  Val  (/)  =  f\e  for  linear  and  exponential 

value  functions,  respectively.  These  distributions  were  chosen  to  provide 
a  reasonable  spread  of  values  for  a  ,  the  decay  rate,  and  /? ,  the  initial 
value  for  each  information  package. 

•  Expected  service  times:  “disparate”  and  “similar.”  For  disparate,  we  use 
the  values  (0.1,  3,  7)  time  units15  as  expected  service  times  for  the 

three  information  package  classes,  respectively,  and  (2,  3,  4)  for 
similar.16 

Additionally,  the  following  factors  are  held  constant: 

•  Number  of  information  package  classes:  3 

•  Number  of  information  packages  generated  per  replication:  9 

•  Number  of  replications  per  run  or  scenario: 17  1000 

Figure  17  through  Figure  21  represent  the  results  of  four  scenarios  (runs)  of  the 
static  simulation — one  for  each  combination  of  value  function  type  (linear,  exponential)18 


15  Minutes  will  be  used  as  the  unit  of  time  in  the  simulations. 

16  These  expected  service  time  values  are  based  in  part  on  the  DARPA  C2  Experiment  data.  See 
Figure  8. 

D  A  run  (or  scenario)  is  defined  by  the  levels  of  the  factors  used.  For  example,  value  function  type  = 
“Exponential”  and  expected  service  times  =  “Similar”  define  Run  1. 
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and  expected  service  times  (similar  or  disparate).  For  each  run,  1000  replications  were 
performed.  For  each  replication,  9  information  packages  of  the  appropriate  type  were 
generated  through  random  draw  of  a  and  /?  from  their  respective  distributions.  One  of 
three  classes  was  (uniformly)  randomly  assigned  to  each  information  package,  along  with 
its  associated  expected  service  time.  During  the  simulation,  an  optimal  processing 
sequence  is  determined  for  each  replication  (optimal  with  respect  to  V'  ),  and  the  value 
of  V'  obtained  by  processing  the  nine  information  packages  in  this  sequence  is  compared 
with  that  obtained  by  processing  the  information  packages  in  sequences  governed  by  each 
of  the  heuristics,  HI  through  H6.19  The  goal  is  to  find  one  or  more  heuristics  that 
perform  in  near-optimal  fashion  with  respect  to  V'  .  The  results  are  presented  below. 


c 

o 

r 

o 

Q. 

O 


□ 

rvhl 

□ 

rvh2 

□ 

< 

=r 

GO 

□ 

rvh4 

□ 

rvh5 

Figure  17 


%  of  Optimal 

Scenario  1A — Exponential  Value  Functions,  Similar  Service  Times 


For  a  particular  run  in  the  comparison  study,  only  a  single  type  of  value  function  is  used 
(exponential  or  linear).  Runs  where  value  function  types  are  “mixed”  are  discussed  in  Section  3  (see  Table 
5). 

I9  The  remaining  heuristics  will  be  examined  in  the  dynamic  simulation. 
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rvh6E 

□  eta=0 

□  eta=.25 

□  eta=.5 

□  eta=.75 

□  eta=1 

'  .r  eta=1 
>T  eta=0 

%  of  Optimal 

Figure  18  Scenario  IB — Exponential  Value  Functions,  Similar  Service  Times 

Figure  17  and  Figure  18  together  show  the  results  from  Scenario  1,  defined  by 
exponential  value  function  type  and  similar  expected  service  times.  These  results  are 
explained  as  follows:  “rvh*”  in  the  chart  legend  stands  for  the  realized  value  histogram 
obtained  by  processing  the  information  packages  in  the  sequence  governed  by  strategy 
FI*.  For  each  heuristic,  a  histogram  is  constructed  that  represents  how  well  the  heuristic 
performed  with  respect  to  the  optimal  sequence.  The  horizontal  axis  represents 
percentage  of  optimal;  the  vertical  axis  represents  the  proportion  of  the  total  number  of 
replications  that  performed  at  the  given  level  (or  “bin”  along  the  horizontal  axis).  For 
example,  Figure  17  shows  a  distribution  for  heuristic  H 1  that  is  centered  at  0.25  with  a 
fairly  large  variance,  while  H5  is  centered  at  a  value  greater  than  0.95  with  a  considerably 
smaller  variance.  Hence,  in  the  vast  majority  of  the  1000  replications,  H5  performed 
nearly  as  well  as  the  exact  optimal  solution.  Heuristic  H6E  (Figure  18)  did  not  perform 
well  for  any  value  of  77 .  It  turns  out  that  H6E  and  H6L  do  not  perform  well  for  any 
combination  of  the  factors  considered,  and,  hence,  will  be  excluded  from  subsequent 
discussion.  Results  for  remaining  runs  are  now  presented  graphically  in  Figure  19 
through  Figure  21.  Results  are  summarized  numerically  in  Table  2. 
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1 


Figure  19  Scenario  2 — Exponential  Value  Functions,  Disparate  Service  Times 


Figure  20 


Scenario  3 — Linear  Value  Functions,  Similar  Service  Times 
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Figure  21  Scenario  4 — Linear  Value  Functions,  Disparate  Service  Times 


Mean  Results  by  Heuristic 

HI 

H2 

H3 

H4 

H5 

Run 

1 

0.263 

0.570 

0.772 

0.934 

0.976 

2 

0.876 

0.305 

0.355 

0.955 

0.986 

3 

0.717 

0.798 

0.821 

0.793 

0.871 

4 

0.886 

0.574 

0.603 

0.787 

0.958 

Table  2  Comparison  Study  Summary  Results 


The  results  in  Table  2  are  computed  in  a  slightly  different  fashion  than  those  in 
Figure  17  through  Figure  21.  The  results  in  these  figures  are  computed  by  normalizing 
each  heuristic  result  against  the  optimal  result  for  each  replication;  this  yields  1000 
observations  for  each  heuristic  for  each  run,  from  which  the  respective  histograms  are 
generated.  In  Table  2,  however,  we  compute  the  mean  in  the  more  traditional  sense  by 
dividing  mean  V'  for  each  heuristic  by  that  for  the  optimal  sequence.  Because  of  the 
large  number  of  replications,  however,  and  due  to  the  nature  of  the  problem,  the  mean 
values  reported  in  Table  2  are  very  close  to  their  respective  empirical  distribution  means 
as  reflected  in  the  histograms  in  Figure  17  through  Figure  21. 
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Based  on  the  results  in  this  limited  examination,  we  conclude  that  H5  performs 
quite  well  when  compared  to  optimal  sequencing  and  other  heuristic  strategies. 
Additionally,  with  the  exception  of  H4,  the  remaining  policies  are  quite  sensitive  to  the 
input  factors,  making  them  less  effective  as  information  package  processing  strategies. 
One  result  of  note,  however,  is  the  effectiveness  of  HI  in  Run  4  (linear  value  functions, 
disparate  service  times);  this  is  to  be  expected,  though,  as  the  conditions  in  Run  4  nearly 
match  the  assumptions  in  Theorem  3.5.  A  more  thorough  examination  is  required  prior 
to  drawing  any  further  conclusions.  This  is  accomplished  in  the  following  sections. 

2.  Validation  of  Comparison  Study 

We  continue  our  investigation  into  the  static  scheduling  problem.  One  issue  that 
must  be  resolved  is  whether  the  results  above  will  hold  in  cases  where  there  are  a  larger 
number  of  information  packages  to  be  sequenced  for  processing.  Here,  we  address  this 
issue  with  a  more  comprehensive  set  of  simulation  runs.  This  set  of  runs  is  described 
below  in  Table  3.  As  in  the  first  set  of  runs,  there  are  two  value  function  types 
considered,  exponential  (E)  and  linear  (L),  and  two  sets  of  expected  service  times — 
similar  (S)  and  disparate  (D).  Numerical  values  for  expected  service  times  are  identical 
to  those  used  previously.  The  number  of  information  packages  (IPs)  alternates  between  9 
and  100,  and  we  allow  two  levels  each  for  a  and  /? — high  (H)  and  low  (L).  As  before, 
each  run  (scenario)  is  replicated  1000  times.  Recall  that  for  exponential  value  functions, 
Val;  (t)  =  fie0'' ,  and  Val  /  (/)  =  /?.  -at  for  linear  value  functions.  Furthermore,  a  and 
/?  are  drawn  from  Gamma  distributions  (unlike  the  previous  set  of  runs)  as  follows: 

•  High  «~T(5,  0.15) 

•  Low  a  ~  T(l,  0.2) 

•  High  ~  E (1.5,  7) 

•  Low  (3  ~  r(i,  4) 
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Run 

Val  Fct 

#IPs 

a 

P 

E[S] 

la 

E 

9 

L 

H 

S 

1b 

E 

100 

L 

H 

s 

2a 

E 

9 

L 

H 

D 

2b 

E 

100 

L 

H 

D 

3a 

E 

9 

L 

L 

S 

3b 

E 

100 

L 

L 

S 

4a 

E 

9 

L 

L 

D 

4b 

E 

100 

L 

L 

D 

5a 

E 

9 

H 

H 

S 

5b 

E 

100 

H 

H 

S 

6a 

E 

9 

H 

H 

D 

6b 

E 

100 

H 

H 

D 

7a 

E 

9 

H 

L 

S 

7b 

E 

100 

H 

L 

S 

8a 

E 

9 

H 

L 

D 

8b 

E 

100 

H 

L 

D 

9a 

L 

9 

L 

H 

S 

9b 

L 

100 

L 

H 

S 

10a 

L 

9 

L 

H 

D 

10b 

L 

100 

L 

H 

D 

11a 

L 

9 

L 

L 

S 

11b 

L 

100 

L 

L 

S 

12a 

L 

9 

L 

L 

D 

12b 

L 

100 

L 

L 

D 

13a 

L 

9 

H 

H 

S 

13b 

L 

100 

H 

H 

S 

14a 

L 

9 

H 

H 

D 

14b 

L 

100 

H 

H 

D 

15a 

L 

9 

H 

L 

S 

15b 

L 

100 

H 

L 

S 

16a 

L 

9 

H 

L 

D 

16b 

L 

100 

H 

L 

D 

Table  3  Validation  Run  Matrix 


Using  these  distributions  for  a  and  /?  allow  for  a  more  realistic  representation  of 
the  types  of  information  packages  present  in  the  physical  system.  For  example,  an 
information  package  with  a  high  a  (value  function  decay  rate)  might  refer  to  a  Time 
Sensitive  Target  (TST);  one  with  a  high  /?  (initial  value)  might  refer  to  a  High  Value 
Target  (HVT).  An  information  package  with  high  a  and  ft  could  refer  to  a  target  that  is 
both  a  TST  and  HVT  (see  [39]  for  target  type  descriptions  and  definitions).  Additionally, 
the  two-parameter  Gamma  distribution  allows  for  tailoring  of  the  “tails”  of  the 
distribution  in  addition  to  the  mean,  as  seen  in  the  high  values  for  a  and  (3  (see  Figure 
22  below).  The  “fatter”  tails  on  the  right  result  in  the  desired  effect  of  a  reasonable 
likelihood  of  very  high  values  (relatively)  for  these  cases. 
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Figure  22  Value  Function  Parameters 


As  before,  an  optimal  solution  is  computed  for  the  runs  having  9  information 
packages,  and  results  for  each  heuristic  are  compared  to  it.  For  runs  involving  100 
information  packages,  exact  optimal  solutions  cannot  be  determined.  H3,  seen  as  a 
“reasonable”  greedy  policy,  is  used  as  a  basis  of  comparison  in  these  cases.  Figure  23 
below  depicts  an  example  of  the  detailed  results  obtained  during  validation  (see 
Appendix  A  for  the  complete  set  of  detailed  results  for  all  16  runs),  and  consolidated 
results  are  shown  in  Table  4. 
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100  IPs 
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rvhl 

□ 

rvh2 
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rvh3 

□ 
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rvh5 

Figure  23  Comparison  Study  Validation  Sample  Results  (Run  2)20 
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-11  Note  that  the  vertical  axes  here  reflect  “count”  (out  of  1000  replications)  rather  than  “proportion.” 
The  interpretation  of  these  graphs  is  the  same  as  before,  however. 
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Heuristic  Mean  Values 

HI 

H2 

H3 

H4 

H5 

Best  9 

Best 

9 

100 

9 

100 

9 

100 

9 

100 

9 

100 

IP 

100  IP 

i 

0.608 

0.119 

0.910 

0.784 

0.926 

1.000 

0.905 

0.978 

0.925 

1.025 

H3 

H5 

2 

0.819 

1.791 

0.724 

0.788 

0.739 

1.000 

0.833 

1.213 

0.972 

2.363 

H5 

H5 

3 

0.585 

0.108 

0.919 

0.810 

0.933 

1.000 

0.916 

0.977 

0.930 

1.022 

H3 

H5 

4 

0.804 

1.654 

0.734 

0.810 

0.749 

1.000 

0.851 

1.220 

0.975 

2.267 

H5 

H5 

5 

0.206 

0.042 

0.657 

0.557 

0.752 

1.000 

0.977 

1.709 

0.989 

1.763 

H5 

H5 

6 

0.859 

22.063 

0.339 

0.675 

0.358 

1.000 

0.949 

25.610 

0.996 

30.720 

H5 

H5 

7 

0.195 

0.037 

0.686 

0.616 

0.763 

1.000 

0.977 

1.676 

0.988 

1.727 

H5 

H5 

c 

3 

8 

0.837 

21.522 

0.371 

0.716 

0.388 

1.000 

0.938 

27.134 

0.996 

32.312 

H5 

H5 

a: 

9 

0.975 

0.595 

0.938 

1.025 

0.929 

1.000 

0.924 

0.992 

0.932 

1.025 

HI 

H5 

10 

0.981 

1.149 

0.894 

1.015 

0.889 

1.000 

0.891 

1.001 

0.961 

1.433 

HI 

H5 

11 

0.863 

0.293 

0.919 

0.992 

0.917 

1.000 

0.904 

0.985 

0.915 

1.021 

H2 

H5 

12 

0.916 

1.092 

0.845 

0.978 

0.846 

1.000 

0.858 

1.007 

0.958 

1.557 

H5 

H5 

13 

0.712 

0.165 

0.918 

0.949 

0.924 

1.000 

0.920 

0.994 

0.941 

1.058 

H5 

H5 

14 

0.854 

1.777 

0.756 

0.940 

0.764 

1.000 

0.812 

1.100 

0.976 

2.338 

H5 

H5 

15 

0.390 

0.069 

0.901 

0.916 

0.919 

1.000 

0.932 

0.994 

0.951 

1.060 

H5 

H5 

16 

0.771 

2.164 

0.649 

0.922 

0.659 

1.000 

0.818 

1.263 

0.985 

2.843 

H5 

H5 

Table  4  Comparison  Study  Validation  Consolidated  Results 


Recall  that  the  exact  optimal  solution  is  the  basis  of  comparison  for  runs 
involving  9  information  packages,  while  H3  serves  this  role  for  runs  having  100 
information  packages  (H3  is  seen  as  a  “reasonable”  greedy  policy,  making  it  a  reasonable 
basis  of  comparison).  Therefore,  the  numerical  results  shown  in  Table  4  can  only  be 
compared  against  one  another  for  cases  where  the  number  of  information  packages  is  the 
same.  While  results  for  the  9-IP  cases  are  strictly  less  than  one  (heuristics  are  inherently 
sub-optimal),  results  for  the  100-IP  cases  can  be  greater  than  one  (indicating  that  the 
heuristic  in  question  performed  better  than  H3)  or  not.  The  columns  labeled  “Best  9  IP” 
and  “Best  100  IP”  indicate  the  best  perfonning  policy  for  the  9-IP  and  100-IP  cases, 
respectively,  for  each  run. 

We  observe  the  following  based  on  these  validation  runs: 


•  If  a  set  of  factors  results  in  H5  outperforming  all  other  heuristics  for  the  9- 
IP  case,  then  H5  outperforms  the  others  in  the  100-IP  case  for  the  same  set 
of  factors  (in  fact,  H5  outperforms  the  others  in  all  100-IP  cases). 

•  The  five  cases  where  H5  did  not  perfonn  the  best  can  be  explained  as 
follows: 
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0  The  factor  levels  for  runs  10  and  11  nearly  match  the  conditions 
required  by  Theorem  3.5,  so  it  is  no  surprise  that  HI  performs  the  best. 

0  For  runs  1,  3,  and  11,  the  differences  in  performance  between  H5  and 
the  highest  performing  heuristic  are  essentially  negligible. 

Additionally,  these  are  all  cases  having  similar  expected  service  times. 
The  physical  system  is  assumed  (with  justification)  to  have 
information  package  classes  with  disparate  expected  service  times 
only.21 

•  The  bottom  line  is  that  there  is  enough  consistency  between  the  9-IP  and 
100-IP  runs  to  warrant  further  investigation  of  H5  as  a  promising 
information  package  processing  assignment  sequencing  heuristic. 

3.  Robustness  of  Policy  H5 

We  now  complete  our  investigation  of  the  static  scheduling  problem  by 
presenting  results  for  cases  where  information  packages  are  no  longer  homogeneous  with 
respect  to  value  function  type;  we  allow  them  to  be  mixed,  as  one  would  expect  to 
encounter  in  the  physical  system.  Additionally,  we  now  include  a  third  value  function 
type,  the  step  function,  defined  as  follows. 


Valq(t)  =  pqH{t)H(rq-t) 

\Pq 

0  otherwise 


(4.2) 


Note  that  H(t)  is  the  unit  (or  Heaviside)  step  function.  This  set  of  simulation  runs 

demonstrates  the  robustness  of  heuristic  H5  to  variations  in  five  significant  factors, 
including  three  used  in  previous  simulation  runs — a  ,  /?,  and  expected  service  time  (all 
defined  as  before) — as  well  as  two  new  factors:  value  function  mix,  which  describes  the 
expected  proportions  of  information  packages  having  exponential  (E),  linear  (L),  and  step 


21  For  example,  typical  information  package  classes  in  the  physical  system  include,  among  others, 
video,  still  imagery,  text,  and  voice.  Clearly,  the  processing  time  distributions  for  each  of  these  classes  will 
differ  significantly. 
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(S)  value  function  types  present  in  the  run,  and  r  ,  which  is  the  time  at  which  step  value 
functions  drop  to  zero.  As  with  a  and  /? ,  we  define  two  levels,  high  and  low,  for  the 
factor  r  as  follows: 

•  High  r  ~  r(2,  20) 

•  Low  r  ~  r(l,  10) 

The  run  matrix  is  shown  in  Table  5  below. 


Run 

Reps 

#  IPs 

ValFctMix 

(E/L/S) 

a 

P 

T 

E[S] 

1 

500 

8 

75/15/10 

L 

L 

H 

D 

2 

500 

8 

75/15/10 

L 

L 

H 

S 

3 

500 

8 

75/15/10 

L 

L 

L 

D 

4 

500 

8 

75/15/10 

L 

L 

L 

S 

5 

500 

8 

75/15/10 

L 

H 

H 

D 

6 

500 

8 

75/15/10 

L 

H 

H 

S 

7 

500 

8 

75/15/10 

L 

H 

L 

D 

8 

500 

8 

75/15/10 

L 

H 

L 

S 

9 

500 

8 

75/15/10 

H 

L 

H 

D 

10 

500 

8 

75/15/10 

H 

L 

H 

S 

11 

500 

8 

75/15/10 

H 

L 

L 

D 

12 

500 

8 

75/15/10 

H 

L 

L 

S 

13 

500 

8 

75/15/10 

H 

H 

H 

D 

14 

500 

8 

75/15/10 

H 

H 

H 

S 

15 

500 

8 

75/15/10 

H 

H 

L 

D 

16 

500 

8 

75/15/10 

H 

H 

L 

S 

17 

500 

8 

33/33/33 

L 

L 

H 

D 

18 

500 

8 

33/33/33 

L 

L 

H 

S 

19 

500 

8 

33/33/33 

L 

L 

L 

D 

20 

500 

8 

33/33/33 

L 

L 

L 

S 

21 

500 

8 

33/33/33 

L 

H 

H 

D 

22 

500 

8 

33/33/33 

L 

H 

H 

S 

23 

500 

8 

33/33/33 

L 

H 

L 

D 

24 

500 

8 

33/33/33 

L 

H 

L 

S 

25 

500 

8 

33/33/33 

H 

L 

H 

D 

26 

500 

8 

33/33/33 

H 

L 

H 

S 

27 

500 

8 

33/33/33 

H 

L 

L 

D 

28 

500 

8 

33/33/33 

H 

L 

L 

S 

29 

500 

8 

33/33/33 

H 

H 

H 

D 

30 

500 

8 

33/33/33 

H 

H 

H 

S 

31 

500 

8 

33/33/33 

H 

H 

L 

D 

32 

500 

8 

33/33/33 

H 

H 

L 

S 

Table  5  Robust  Static  Scheduling  Simulation  Run  Matrix 
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Note  that  [J  represents  the  initial  value  for  all  three  types  of  value  functions,  and 
a  represents  the  decay  rate  for  linear  and  exponential  value  functions.  Note  also  that, 
unlike  previous  runs,  500  replications  are  executed  per  run,  with  8  information  packages 
per  replication.  This  decrease  (from  1000  and  9,  respectively)  is  due  to  run  time 
considerations  only,  and  does  not  significantly  affect  the  results.  As  in  previous  runs,  the 
result  for  each  heuristic  is  compared  to  the  exact  optimal  solution.  The  results  for  these 
runs  are  shown  in  Table  6.  Each  result  represents  an  average  over  the  500  replications 
for  each  run. 


Run 

H2 

H3 

H4 

H5 

Spread 

1 

0.815 

0.818 

0.861 

0.964 

0.150 

2 

0.924 

0.924 

0.892 

0.907 

0.032 

3 

0.783 

0.792 

0.853 

0.967 

0.183  i 

4 

0.909 

0.921 

0.907 

0.917 

0.014 

5 

0.802 

0.809 

0.840 

0.955 

0.153 

6 

0.915 

0.915 

0.882 

0.894 

0.032 

7 

0.754 

0.765 

0.819 

0.958 

0.204  i 

8 

0.891 

0.903 

0.889 

0.908 

0.019  I 

9 

0.652 

0.669 

0.872 

0.975 

0.323 

10 

0.835 

0.880 

0.894 

0.910 

0.074 

11 

0.577 

0.603 

0.882 

0.980 

0.403 

12 

0.752 

0.837 

0.916 

0.929 

0.177  i 

13 

0.629 

0.657 

0.827 

0.971 

0.342 

14 

0.828 

0.882 

0.878 

0.897 

0.069 

15 

0.575 

0.606 

0.861 

0.976 

0.400 

16 

0.765 

0.857 

0.901 

0.913 

0.148 

17 

0.863 

0.863 

0.871 

0.952 

0.090  : 

18 

0.910 

0.909 

0.892 

0.896 

0.018 

19 

0.791 

0.806 

0.855 

0.956 

0.165 

20 

0.855 

0.873 

0.890 

0.891 

0.036 

21 

0.872 

0.870 

0.874 

0.949 

0.079  i 

22 

0.907 

0.902 

0.883 

0.893 

0.024 

23 

0.791 

0.802 

0.848 

0.944 

0.153 

24 

0.847 

0.862 

0.869 

0.882 

0.035 

25 

0.774 

0.784 

0.840 

0.953 

0.178 

26 

0.883 

0.897 

0.875 

0.886 

0.021 

27 

0.651 

0.689 

0.833 

0.961 

0.310 

28 

0.774 

0.838 

0.880 

0.887 

0.112 

29 

0.769 

0.786 

0.819 

0.945 

0.176 

30 

0.868 

0.891 

0.882 

0.892 

0.024  | 

31 

0.655 

0.695 

0.808 

0.955 

0.300  i 

32 

0.773 

0.841 

0.879 

0.886 

0.114 

mean 

0.793 

0.817 

0.868 

0.930 

0.142 

std  dev 

0.099 

0.092 

0.027 

0.033 

0.117  | 

Table  6  Robust  Static  Simulation  Results — Average  Percentage  of  Optimal 
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Note  that  HI  is  not  included  in  the  above  results;  it  is  not  defined  for  step 
functions  (there  is  no  decay  rate  for  step  functions).  The  values  in  this  table  represent  the 
mean  realized  value  obtained  (over  the  500  replications)  for  each  policy  for  each  run, 
expressed  as  a  proportion  of  the  optimal  result  (for  each  run).  The  column  titled 
“Spread”  records  the  difference  between  the  largest  and  smallest  mean  result  for  each 
run.  These  results  indicate  that  H5  is  quite  robust  to  variations  in  the  significant  factors. 
On  average,  sequencing  information  packages  according  to  H5  results  in  a  V'  that  is 
only  7%  less  than  that  obtained  from  optimal  sequencing.  The  six  runs  for  which  H5 
does  not  perform  best  all  have  very  small  spreads,  indicating  that  the  four  heuristics  are 
nearly  indistinguishable  in  terms  of  their  perfonnance.  These  results  clearly  suggest  that 
H5  is  robust  to  variations  in  the  significant  factors. 

D.  DYNAMIC  SCHEDULING  SIMULATION 

The  dynamic  scheduling  simulation  is  designed  to  be  a  credible  emulation  of  the 
physical  system,  with  random  arrival  and  service  times  and  a  realistic  mix  of  value 
function  types.  Its  purpose  is  to  compare  promising  heuristic  sequencing  strategies  to 
each  other  and,  ultimately,  to  a  representation  of  the  “status  quo” — how  these  processing 
decisions  are  currently  made.  Unlike  the  static  scheduling  simulation,  which  represents 
only  a  snapshot  in  time,  the  dynamic  simulation  represents  the  entire  duration  of  an 
operation.  We  introduce  a  new  heuristic  called  “Status  Quo”  as  an  estimate  of  how 
information  package  assignment  decisions  are  currently22  made.  Status  Quo  is  defined  to 
be  a  weighted  combination  of  three  previously  defined  heuristics  as  follows:  60%  LIFO, 
30%  H3,  and  10%  Random  [40].  Status  Quo  is  executed  in  the  simulation  by  randomly 
choosing  one  of  its  three  component  heuristics  (according  to  their  respective  weights)  to 
use  each  time  an  assignment  decision  is  made.  The  goal  is  to  measure  how  the  policies 
previously  defined  (specifically  H5)  compare  with  Status  Quo  using  Total  Realized 
Value  as  the  metric. 


22  As  observed  in  the  2006  DARPA  Command  and  Control  experiments. 
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The  dynamic  simulation  starts  with  zero  information  packages  in  Oval  2. 
Information  packages  begin  arriving  according  to  a  process  assumed  to  be  Poisson  with 
rate  i  =  1,2,..., NumClass  and  are  randomly  assigned  a  value  function  type  and 

associated  parameters  (a ,  fJ ,  and/or  z  ).  Information  packages  are  assigned,  without 
preemption,  to  available  processors  as  per  the  method  described  in  Figure  16,  with 
service  times  assumed  to  follow  a  Poisson  process  with  rate  i  =  1,2,..., NumClass  ,23 

This  process  is  repeated  in  parallel,  using  identical  sets  of  information  packages,  for  each 
heuristic  (H3,  H4,  H5,  LIFO,  FIFO,  Random,  and  StatusQuo).  The  simulation  ends  when 
the  last  information  package  arrives,  leaving  many  infonnation  packages  still 
unprocessed — a  realistic  outcome  mirroring  the  physical  system.  Realized  value  is 
computed  for  each  information  package  processed  using  the  expression 

VJr=Valj(zj+sj),  (4.3) 

and  summed  to  detennine  Total  Realized  Value.  This  value  is  recorded  for  each 
heuristic.  MATLAB  code  for  the  dynamic  simulation  can  be  found  in  Appendix  A. 
Figure  24  shows  a  typical  outcome  of  a  single  replication  of  the  dynamic  simulation, 
revealing  how  realized  value  grows  over  the  duration  of  the  simulation  (operation) 
according  to  the  processing  strategy  utilized. 


23  The  Poisson  assumptions  are  reasonable  based  on  data  obtained  from  the  DARPA  Command  and 
Control  experiments.  See  Figure  7  and  Figure  8. 
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Value 


Realized  Value  vs.  Time 


Figure  24  Typical  Dynamic  Simulation  Replication  Results 

A  full-factorial  design  is  executed  using  the  following  seven  factors,  each  having 
two  levels  (for  a  total  of  27  =  128  runs): 

•  Number  of  servers  (processors) — 3  or  5. 

•  Value  Function  Mix — 33/33/33  or  75/15/10  for  exponential,  linear,  and 
step  value  function  types,  respectively. 

•  a  (linear  and  exponential  value  functions) — low  (L)  and  high  (FI),  as  in 
the  static  scheduling  simulation  runs. 

•  / 3  (linear  and  exponential  value  functions) — low  (L)  and  high  (FI),  as  in 
the  static  scheduling  simulation  runs. 

•  / 3  (step  functions) — unlike  the  static  case,  we  allow  the  initial  value  for 
step  functions  to  differ  from  that  for  linear  and  exponential  value 
functions.  Two  levels  are  used  as  follows. 

0  High  4) 

0  Low  /?~T(0.5,  2) 
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•  t  (step  functions) — low  (L)  and  high  (H),  as  in  the  static  simulation  runs. 

•  A  (arrival  rate) — two  sets  of  arrival  rates  are  used:  [0.631,  0.5,  0.8]24  and 
[1.2,  0.2,  0.6]  arrivals  per  minute  for  infonnation  packages  of  class  1,  2, 
and  3,  respectively. 

In  addition  to  these  varying  factors,  the  following  factors  are  held  constant  for 
each  run. 

•  Number  of  replications — 250 

•  Number  of  information  packages — 1000 

•  Number  of  information  package  classes — 3 

•  /j  (processing  rates) — one  set  of  processing  rates  is  used:  [0.134,  0.1,  l]25 
processing  completions  per  minute  for  information  packages  of  class  1,  2, 
and  3,  respectively. 

For  each  of  the  128  runs,  Total  Realized  Value  is  computed  for  each  policy  for 
each  replication,  and  then  averaged  over  all  replications.  The  detailed  run  matrix  can  be 
found  in  Appendix  A,  and  the  results  are  summarized  in  Table  7  below. 


24  This  set  of  arrival  rates  is  based  in  part  on  DARPA  C2  Experiment  data. 
2^  Service  rates  are  also  based  in  part  on  DARPA  C2  Experiment  data 
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Mean  Values  for  each  heuristic  for  each  run 

Run 

H5  |  H4  |  H3 

LIFO 

Random 

StatusQuo 

FIFO 

h5/StatusQuo 

1 

1556 

1424 

1357 

608 

166 

887 

92 

1.75 

2 

1434 

1334 

1280 

533 

151 

815 

82 

1.76 

3 

1629 

1476 

1407 

632 

186 

915 

98 

1.78 

4 

1493 

1368 

1324 

557 

164 

837 

85 

1.78 

5 

1935 

1750 

1622 

781 

192 

1013 

98 

1.91 

6 

1765 

1616 

1523 

696 

174 

932 

91 

1.89 

7 

2288 

2009 

1923 

894 

274 

1245 

131 

1.84 

8 

2093 

1874 

1807 

810 

239 

1147 

117 

1.82 

9 

4094 

3691 

3553 

1567 

561 

2368 

340 

1.73 

10 

3798 

3481 

3385 

1395 

512 

2241 

312 

1.69 

11 

4161 

3716 

3567 

1605 

582 

2411 

354 

1.73 

12 

3847 

3492 

3403 

1439 

527 

2235 

316 

1.72 

13 

4415 

3879 

3657 

1751 

594 

2490 

355 

1.77 

14 

4028 

3610 

3449 

1553 

530 

2310 

313 

1.74 

15 

4673 

4014 

3821 

1876 

672 

2629 

389 

1.78 

16 

4253 

3744 

3601 

1686 

606 

2457 

345 

1.73 

17 

821 

795 

730 

357 

37 

372 

15 

2.21 

18 

714 

693 

656 

304 

31 

336 

13 

2.12 

19 

907 

852 

806 

384 

57 

438 

22 

2.07 

20 

823 

782 

757 

333 

51 

386 

21 

2.13 

21 

1232 

1183 

1090 

541 

63 

557 

25 

2.21 

22 

1099 

1051 

1014 

460 

54 

499 

22 

2.20 

23 

1664 

1540 

1482 

652 

137 

878 

55 

1.89 

24 

1547 

1442 

1396 

576 

120 

810 

47 

1.91 

25 

2512 

2379 

2163 

1022 

144 

1276 

57 

1.97 

26 

2267 

2160 

2032 

874 

125 

1160 

43 

1.95 

27 

2580 

2405 

2215 

1028 

165 

1283 

66 

2.01 

28 

2387 

2227 

2111 

912 

146 

1214 

58 

1.97 

29 

2885 

2700 

2418 

1201 

171 

1405 

66 

2.05 

30 

2612 

2446 

2260 

1046 

146 

1309 

54 

1.99 

31 

3224 

2908 

2692 

1323 

236 

1587 

90 

2.03 

32 

2970 

2714 

2551 

1162 

214 

1505 

82 

1.97 

33 

1648 

1477 

1362 

680 

165 

896 

88 

1.84 

34 

1486 

1348 

1279 

595 

145 

802 

75 

1.85 

35 

1640 

1468 

1349 

665 

159 

878 

85 

1.87 

36 

1487 

1355 

1279 

592 

143 

811 

76 

1.83 

37 

1853 

1641 

1528 

757 

196 

991 

99 

1.87 

38 

1691 

1523 

1443 

669 

170 

925 

89 

1.83 

39 

1745 

1560 

1420 

728 

171 

921 

89 

1.89 

40 

1577 

1435 

1337 

630 

147 

850 

76 

1.85 

41 

4286 

3729 

3425 

1773 

499 

2286 

272 

1.88 

42 

3824 

3391 

3189 

1557 

432 

2086 

242 

1.83 

43 

4237 

3715 

3403 

1774 

472 

2282 

266 

1.86 

44 

3830 

3403 

3193 

1564 

424 

2109 

240 

1.82 

45 

4427 

3830 

3510 

1858 

509 

2352 

284 

1.88 

46 

3958 

3479 

3245 

1653 

453 

2161 

248 

1.83 

47 

4354 

3774 

3448 

1814 

481 

2336 

269 

1.86 

48 

3878 

3442 

3192 

1597 

429 

2127 

238 

1.82 

49 

692 

664 

557 

312 

31 

311 

13 

2.22 

50 

597 

571 

512 

275 

27 

279 

12 

2.14 
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Mean  Values  for  each  heuristic  for  each  run 

Run 

H5  |  H4  |  H3 

UFO 

Random 

StatusQuo 

FIFO 

h5/StatusQuo 

51 

673 

656 

536 

317 

24 

300 

10 

2.24 

52 

569 

556 

481 

266 

22 

268 

10 

2.12 

53 

943 

896 

806 

402 

55 

490 

22 

1.93 

54 

832 

792 

741 

348 

47 

439 

21 

1.90 

55 

813 

788 

673 

372 

35 

365 

16 

2.23 

56 

684 

662 

600 

309 

27 

317 

11 

2.16 

57 

1980 

1882 

1547 

863 

92 

954 

37 

2.08 

58 

1697 

1622 

1408 

741 

75 

831 

34 

2.04 

59 

1985 

1869 

1573 

879 

98 

925 

43 

2.15 

60 

1686 

1608 

1409 

740 

82 

852 

37 

1.98 

61 

2081 

1977 

1624 

926 

102 

994 

45 

2.09 

62 

1792 

1701 

1479 

781 

83 

896 

34 

2.00 

63 

2202 

2056 

1774 

969 

123 

1083 

52 

2.03 

64 

1927 

1791 

1614 

825 

102 

968 

42 

1.99 

65 

1757 

1702 

1662 

983 

312 

1243 

183 

1.41 

66 

1679 

1631 

1615 

880 

276 

1180 

167 

1.42 

67 

1866 

1796 

1768 

1032 

347 

1318 

198 

1.42 

68 

1764 

1706 

1689 

928 

310 

1240 

178 

1.42 

69 

2217 

2136 

2093 

1282 

369 

1519 

207 

1.46 

70 

2071 

2000 

1980 

1146 

316 

1402 

177 

1.48 

71 

2644 

2510 

2463 

1476 

515 

1848 

276 

1.43 

72 

2503 

2389 

2358 

1332 

453 

1720 

236 

1.46 

73 

4704 

4518 

4443 

2582 

1047 

3414 

701 

1.38 

74 

4451 

4291 

4246 

2281 

909 

3173 

595 

1.40 

75 

4749 

4532 

4454 

2618 

1052 

3428 

692 

1.39 

76 

4491 

4315 

4279 

2346 

958 

3236 

615 

1.39 

77 

5076 

4821 

4657 

2883 

1104 

3608 

714 

1.41 

78 

4808 

4583 

4471 

2581 

969 

3373 

618 

1.43 

79 

5412 

5075 

4940 

3056 

1217 

3797 

763 

1.43 

80 

5084 

4798 

4712 

2768 

1116 

3572 

696 

1.42 

81 

937 

926 

917 

571 

76 

577 

31 

1.62 

82 

839 

829 

844 

488 

63 

513 

25 

1.64 

83 

1063 

1037 

1039 

622 

113 

683 

50 

1.56 

84 

956 

937 

955 

542 

94 

609 

41 

1.57 

85 

1905 

1839 

1853 

1081 

279 

1319 

124 

1.44 

86 

1778 

1721 

1745 

937 

234 

1214 

98 

1.46 

87 

1409 

1384 

1405 

883 

125 

874 

54 

1.61 

88 

1280 

1259 

1313 

765 

102 

808 

40 

1.58 

89 

2972 

2884 

2804 

1691 

331 

1950 

141 

1.52 

90 

2745 

2665 

2642 

1501 

276 

1808 

119 

1.52 

91 

2885 

2832 

2718 

1653 

289 

1906 

122 

1.51 

92 

2673 

2614 

2551 

1461 

244 

1755 

102 

1.52 

93 

3733 

3564 

3484 

2118 

488 

2477 

210 

1.51 

94 

3464 

3322 

3258 

1897 

419 

2278 

182 

1.52 

95 

3299 

3200 

3055 

1925 

343 

2140 

148 

1.54 

96 

3036 

2959 

2886 

1711 

293 

1991 

118 

1.52 

97 

1923 

1836 

1761 

1090 

312 

1290 

177 

1.49 

98 

1782 

1712 

1669 

979 

276 

1206 

158 

1.48 

99 

1894 

1816 

1737 

1087 

303 

1272 

174 

1.49 

100 

1765 

1698 

1652 

968 

268 

1181 

152 

1.49 

73 


Mean  Values  for  each  heuristic  for  each  run 

Run 

H5 

H4 

H3 

UFO 

Random 

StatusQuo 

FIFO 

h5/StatusQuo 

101 

2159 

2052 

1982 

1235 

364 

1466 

201 

1.47 

102 

2012 

1925 

1878 

1109 

319 

1372 

175 

1.47 

103 

2037 

1951 

1871 

1178 

317 

1349 

177 

1.51 

104 

1873 

1797 

1751 

1045 

279 

1253 

154 

1.49 

105 

4984 

4699 

4494 

2886 

909 

3376 

554 

1.48 

106 

4645 

4411 

4268 

2583 

809 

3132 

488 

1.48 

107 

4944 

4687 

4471 

2891 

900 

3355 

549 

1.47 

108 

4617 

4392 

4235 

2547 

792 

3123 

473 

1.48 

109 

5058 

4778 

4532 

2951 

911 

3418 

551 

1.48 

110 

4681 

4437 

4274 

2628 

799 

3168 

478 

1.48 

111 

5170 

4861 

4622 

3021 

947 

3510 

569 

1.47 

112 

4776 

4506 

4340 

2675 

843 

3236 

497 

1.48 

113 

811 

801 

729 

502 

56 

474 

27 

1.71 

114 

694 

691 

643 

421 

42 

409 

19 

1.70 

115 

842 

829 

769 

522 

63 

506 

29 

1.66 

116 

729 

716 

683 

437 

52 

443 

25 

1.65 

117 

951 

940 

879 

592 

67 

581 

30 

1.64 

118 

827 

817 

793 

506 

54 

501 

24 

1.65 

119 

1092 

1066 

1012 

645 

111 

720 

51 

1.52 

120 

985 

960 

936 

563 

95 

660 

42 

1.49 

121 

2327 

2280 

2043 

1395 

179 

1442 

82 

1.61 

122 

2035 

1991 

1839 

1186 

148 

1278 

66 

1.59 

123 

2337 

2286 

2073 

1422 

196 

1437 

89 

1.63 

124 

2044 

2003 

1856 

1205 

159 

1255 

71 

1.63 

125 

2433 

2380 

2154 

1497 

194 

1531 

86 

1.59 

126 

2160 

2114 

1968 

1281 

168 

1343 

71 

1.61 

127 

2582 

2495 

2299 

1544 

249 

1642 

119 

1.57 

128 

2275 

2212 

2088 

1321 

205 

1437 

88 

1.58 

Table  7  Dynamic  Simulation  Results  Summary 


One  notable  observation  is  that  H5  results  in  the  highest  average  total  realized 
value  for  126  out  of  the  128  runs  (note  green  shading  indicates  the  highest  value  for  each 
run).  As  before,  in  the  cases  where  H5  is  not  the  highest  perfonning  heuristic,  it  is 
virtually  indistinguishable  from  the  one  that  is  (H3  in  both  cases).  H5  also  significantly 
outperforms  Status  Quo  on  every  run — by  as  small  a  factor  as  1.38  (run  73)  and  as  large  a 
factor  as  2.24  (run  51).  Additionally,  we  note  that  H5  performs  well  regardless  of  arrival 
rate,  as  long  as  the  aggregate  arrival  rate  is  much  larger  than  aggregate  service  rate. 
Finally,  implementing  H5  does  not  require  knowledge  of  the  service  time  distributions — 
all  that  is  required  is  mean  service  times.  We  conclude  that  H5  is  a  viable  heuristic  to  use 
in  the  physical  system. 
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E.  IMPLEMENTATION  OF  HEURISTIC  STRATEGY  H5 

Previously  in  this  chapter,  we  have  demonstrated  quite  conclusively  that  H5  is  a 
well-performing  heuristic  strategy  to  employ  when  making  processor  assignment 
decisions  within  the  context  of  the  physical  system.  In  this  section,  we  describe  a  lower 

bound  based  on  H5  and  depict  an  0[n2  j  algorithm  for  implementing  H5. 

1.  Using  H5  to  Construct  a  Lower  Bound  on  V  * 

In  Chapter  III,  we  alluded  to  a  lower  bound  for  optimal  realized  value  V* 
constructed  using  the  function  T.  (t)  .  Here,  we  combine  these  functions  with  strategy  H5 
to  construct  a  lower  bound. 

Let  si  be  the  deterministic  service  time  for  information  package  i.  Define 
Val (t) 

T.  ( t )  = - as  before,  where  Valt  (/ )  is  of  arbitrary  (or  “mixed”)  form,  and  define  the 

Si 

permutation  mapping  cpH5  (y)  as  the  infonnation  package  to  be  processed  /h,  where 
processing  sequence  is  determined  by  heuristic  strategy  H5.  Let 

tf5  =  s  +5  H - \-s  Define  VH5 — the  total  realized  value  obtained  by 

'  <pH\\)  <pH5(2)  <p  (Q)  j 

processing  the  Q  information  packages  according  to  heuristic  strategy  H5 — as  follows: 

r"5=£vwW'f)  <44) 

M 

We  employ  the  same  geometric  reasoning  in  constructing  (4.4)  as  we  did  in  (3.45). 
Finally,  it  is  clear  that  no  heuristic  strategy  can  outperform  an  optimal  strategy,  which 
leads  us  to  the  following  proposition  regarding  a  lower  bound  on  V  * . 
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Proposition  4.1:  VH5  <  V  *  . 

Given  the  empirical  evidence  of  the  perfonnance  of  H5  with  respect  to  optimal 
sequencing,  we  conjecture  that  VH5  is  a  good  lower  bound  for  V  * . 


An  O 


( - 


(n  +  1) 


Implementation  Algorithm  for  H5 


In  this  section,  we  develop  an  O 


i(n  + 1) 


algorithm  for  implementing  H5  for 


the  case  where  all  value  functions  are  exponential,  service  times  are  detenninistic,  and 
m  =  1  server. 


To  begin,  we  define  the  term  f  as  follows: 


_Vali(x,+si) 


(4.5) 


where  Valt  (/)  is  the  value  function  associated  with  information  package  i,  i  =  1,2,..., n  , 
Xj  is  the  time  information  package  i  is  assigned  to  the  processor,  and  si  is  the  service 
time  for  information  package  i  (assumed  detenninistic).  In  Section  B  of  this  chapter,  we 

ValI(E[Di\xl=td~]) 

defined  the  heuristic  H5  based  on  the  sequencing  criterion  - - - - — = - .  An 

£[S,] 

VallE[x,+S,]) 

equivalent  representation  of  this  sequencing  criterion  is  - - — ; ^ - -,  which,  for 


deterministic  service  times  is  equivalent  to  Tj  in  (4.5). 


~  .  Val  (t  +  5.) 

Let  7^  (t)  = - — - —  be  a  continuous  function  of  time  that  represents  the  value 


of  the  H5  sequencing  criterion  for  information  package  i  for  any  assignment  time  t. 
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T.  (?)  is  related  to  the  function  Tt  (?)  introduced  earlier  (Chapter  III,  Section  D,  and  again 
in  Section  El  above);  specifically, 


f,(/)  =  T(l  +  s,). 


(4.6) 


Expression  (4.6)  is  used  to  assign  information  packages  to  processors  as  follows:  at  time 
td  (the  jth  decision  time)  the  information  package  i  that  satisfies  (4.7)  (and  has  not 

already  been  sent  to  a  processor)  is  sent  to  the  available  processor.  This  is  equivalent  to 
assigning  information  packages  according  to  strategy  H5. 


max|f ('•<■)} 


(4.7) 


We  now  construct  a  geometric  argument  which  becomes  the  basis  for  our 
implementation  algorithm.  Figure  25  shows  some  example  functions, 
Ti  (?) ,  i  =  1, 2, . . . ,  4 .  We  begin  by  indexing  the  information  packages  in  decreasing  order 

of  T  (0) .  Recalling  that  T  (?)  is  the  H5  sequencing  criterion,  we  note  that,  assuming 
tdi  =  0 ,  the  first  information  package  chosen  for  processing  is  infonnation  package  1  (i.e., 
^(1)  =  1 )  since  7J  (0)  is  greater  than  the  remaining  initial  values.26  To  determine  <p(  2) 

(i.e.,  the  information  package  that  will  be  processed  second),  we  consider  information 
package  2  as  the  first  candidate,  since  it  has  the  next  highest  initial  value.  We  define  ?( k 

as  the  time  when  f.  (?)  and  fk  (?)  intersect;  see  Figure  25.  We  must  check  the 

intersection  times  of  the  function  T  (?)  for  the  candidate  information  package  with  that 

of  all  remaining  information  packages,  and  compare  these  times  with  td  .  For  our 


26  In  fact,  the  way  we  have  chosen  to  sequence  information  packages  will  always  lead  to  package 
number  one  being  processed  first. 
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example,  we  compute  t23  and  t24.  If  tdi  <  min|?23,t24| ,  then  <p[ 2)  =  2;  otherwise,  we 
re-designate  information  package  3  as  the  candidate  for  (p[l) ,  and  repeat  this  process  of 
comparing  intersection  times  with  td  until  <p(2)  is  determined.  We  formalize  this 
notion — that  of  using  intersection  points  of  the  functions  Tf  (t)  to  determine  processing 
sequence — in  the  algorithm  below. 


Figure  25  H5  Implementation  Algorithm — a  Geometric  Perspective 


First,  we  define  some  terms.  Let  Valt  (/  )  =  f^e  ~a,t  be  the  value  function  for 
information  package  i,  i  =  1, 2, . . . ,  n  ,  let  m  =  1  processor,  and  let  sf  be  the  deterministic 

~  .  Val  it  +  s, )  .  . 

service  time  for  information  package  i.  Define  F  (/)  = - - - - ,  and  let  <p\  j)  be  the 

Si 

information  package  to  be  processed /h. 

We  now  define  the  H5  implementation  algorithm  as  follows. 
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•  Step  1:  Initialize  and  determine  first  information  package.  This  step 
requires  n  function  evaluations. 

0  for  /  =  1  to  n,  evaluate  Ti  (0) ,  end. 

0  (p{\)  =  i :  ft  (0)  >  fk  (0),  \/k  ^  i  (the  first  information  package 

assigned). 


•  Step  2:  for  /  =  2  to  n,  do: 


0  “Discard”  information  package  (p{j  - 1).  (note:  now  n-(j-l) 
information  packages  remain) 

0  Let  c(j,l)  =  i:T.  (0)  >  fk  (0),  \/k  ^  i,  \/i,k  <£  j^(l),...,^(y  -l)}  .  This 

is  the  first  candidate  for  the  /h  assignment;  it  has  the  highest  initial 
value  of  those  information  packages  not  yet  assigned. 

j- 1 

0  Let  td.  =’^s<p(k)  (the/11  decision  time) 

'  k= 1 


0 


Let  t.k  = 


_  \Pthj 


In 


+  aisi-aksk 


ctk~ai 


functions  7) (/)  and  Tk(t )) 


(the  time  of  intersection  of  the 


0  Compare  td  to  tc^  ]]k  for  all  k^c[j,  1)  (each  comparison  requires 
one  function  evaluation) 

if  tdJ  <  UjMk  yk  *  CU> 0 »  then  <pU)  =  cU> !) 

else,  c(j,2)  =  k:tc{jl)k<tdj 

-  repeat  (up  to  n-j  + 1  times)  until  <p(j)  is  determined. 


Step  1  requires  n  function  evaluations  (to  find  the  n  initial  values)  and  results  in 
the  detennination  of  (p{  \  ) .  The  first  loop  through  Step  2  (that  is,  the  loop  during  which 

(p{  2  )  is  determined)  requires  the  examination  of  at  most  n  - 1  intersections,  with  each 

intersection  examination  requiring  one  function  evaluation.  The  second  loop  through 
Step  2  requires  the  examination  of  at  most  n  -  2  intersections  (and  the  execution  of  n  —  2 
function  evaluations),  and  so  on.  So,  completely  sequencing  all  n  infonnation  packages 
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/  \  /  \  n(n  + 1) 

requires  at  most  n  +  (n-\)  +  (n-2)-\ - h 2  + 1  =  — - - -  function  evaluations.  Hence, 


we  conclude  that  this  algorithm  is  O 
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V.  SUMMARY  AND  CONCLUSION 


A.  RESULTS  SUMMARY 

The  results  and  contributions  of  this  research  can  be  categorized  into  three  main 
areas:  1)  development  of  mathematical  concepts  and  tools  that  can  be  used  to  model  the 
dynamics  of  battlefield  information  flow;  2)  analysis  of  these  models  to  derive  optimal 
control  strategies  and  complexity  results;  and  3)  development  of  robust  heuristic  control 
strategies  and  the  construction  of  simulation  platforms  for  strategy  comparison.  We  now 
present  a  detailed  summary  of  these  contributions  and  results. 

1.  Modeling 

•  Some  fundamental  concepts  are  introduced,  namely  the  incorporation  of 
information  value  and  infonnation  volume  (workload)  in  a  dynamic 
information  flow  context. 

•  A  controlled  queue  model  is  developed  for  maximizing  cumulative 
expected  realized  value.  This  involved: 

0  Applying  GIGIm  multi-class  queuing  system  as  the  framework  to 
describe  the  system. 

0  Defining  various  types  of  value  functions  for  information  packages. 

0  Developing  the  optimal  assignment  problem  for  the  queuing  system. 

•  A  discrete  time  information  flow  model  is  developed  to  represent  the 
dynamics  of  total  volume  (or  workload)  in  the  system. 

2.  Analysis  and  Complexity  Results 

•  The  following  results  are  proved  regarding  the  complexity  of  special  cases 
of  the  controlled  queue  model. 
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0  Step 
Val  (t)  = 


Value  Functions  (“single 

\Pq  aq^t^aq^Tq 

I  0  otherwise 


step”):  if 

then  the  optimal  control  problem 


max  E  [V 1  =  max  V  Val  t,  +  V 


x 

/?=i  y 


is  NP-hard. 


0  The  same  is  true  regarding  “multi-step”  value  functions. 


•  The  following  analytical  solutions  for  special  cases  of  the  controlled  queue 
model  are  proved. 

0  Linear  value  functions  ( Valq  (  t  )  =  /3q  -  aqt ) 

-  Proved  that  the  optimal  sequencing  problem 


max 


max  E 


?=i  V 


J 


(5.1) 


is  equivalent  to  the  Quadratic  Assignment  Problem 

Q  Q  Q  Q 

minimize 

1=1  7=1  k= 1  /= 1 
0 

subject  to  Z-fy  =1’ 

i= 1 

Z y<t=l>  v/ 

y'=i 

Ty  e  {0,1} ,  Vz',7 


(5.2) 


The  quadratic  assignment  problem  in  general  is  NP- 
complete,  but  due  to  the  special  structure  of  problems  (5.1) 
and  (5.2),  an  analytic  solution  is  proved.  The  solution  to 
both  is  to  process  information  packages  in  decreasing  order 

of  «./*[>.]• 

This  result  reveals  a  previously-undiscovered  case  of  the 
quadratic  assignment  problem  that  has  an  analytic  optimal 
solution. 

0  Exponential  Value  Functions — the  following  results  are  proved: 
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If  Valq  {t)  =  pqe  at  and  Sq  ~  Exponential  (//J  then 

processing  information  packages  in  decreasing  order  of 
fiq/jq  is  optimal. 

If  Valq  {t)  =  Pqe  a,,t  and  Sq  ~  Exponential  (//)  then 
processing  information  packages  in  decreasing  order  of 

Pq<*q  ■  , 

-  .,  is  optimal. 

(«,+/*) 


If  Valq  (  t  )  =  Pqe  u,,t  and  Sq  ~  Uniform  (a,  h)  then 
processing  information  packages  in  decreasing  order  of 


Pa 


f  -a„a 

e  q  -e 


a„ 


is  optimal.. 


v  ?  J 

Upper  and  lower  bounds  for  the  optimal  solution  in  the  mixed  value 
function  case  are  derived. 


3.  Simulation 


•  Due  to  the  complexity  of  the  problem,  a  set  of  heuristic  information 

package  processing  strategies  is  developed. 

•  Static  scheduling  simulation 

0  Constructed  a  Monte  Carlo  simulation  in  Matlab  to  compare  heuristic 
policies  to  each  other  (in  all  cases)  and  to  exact  optimal  solutions  (for 
small  cases). 

0  Developed  a  heuristic  policy  that  is  robust  to  variations  in  value 
function  type,  value  decay  rate,  initial  value,  and  expected  service 
times. 

•  Dynamic  scheduling  simulation. 

0  Constructed  a  discrete  event  simulation  in  Matlab  to  represent  the 
dynamics  of  the  physical  system. 

0  Developed  a  mathematical  representation  of  “status  quo”  in  battle  field 
information  flow  control.  Integrated  this  model  into  the  simulation 
framework  for  comparison  purposes. 
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0  Developed  a  heuristic  policy  that  significantly  outperfonns  the  status 
quo  and  is  robust  to  variations  in  value  function  type,  value  decay  rate, 
initial  value,  and  arrival  rates. 

0  Developed  an  0(  n2  ]j  algorithm  for  implementing  this  robust  heuristic 
policy. 

These  results  constitute  an  important  initial  step  toward  greater  understanding  of 
the  overall  problem  of  battlefield  information  flow  through  mathematical  modeling  and 
analysis.  More  is  required,  however,  before  we  have  a  representation  of  the  physical 
system  that  is  realistic  and  accurate  enough  for  real-world  application.  In  the  following 
section,  we  describe  some  potential  future  directions  for  this  research. 

B.  FUTURE  WORK 

1.  “Ensembling”  of  Information  Packages 

One  of  the  assumptions  we  have  made  is  that  information  package  value  is 
linearly  additive — meaning  that  the  value  gained  by  processing  a  particular  information 
package  is  independent  of  the  information  packages  that  have  already  been  processed. 
We  realize  that  this  is  not  the  case  in  the  physical  system.  In  fact,  the  value  of  a 
particular  information  package  is  dependent  upon  what  has  already  been  processed  (and 
even,  to  some  extent,  upon  what  is  planned  to  be  processed).  It  is  likely  the  case  that 
several  related  information  packages  (related  in  the  sense  that  they  describe  the  same 
target  or  set  of  targets,  or  the  same  geographical  region)  exist  in  Oval  2 — collected  by 
multiple  sources — and  that  the  value  of  processing  them  in  aggregate  is  greater  than  the 
sum  of  their  individual  values.  This  is  analogous  to  the  advantage  of  multi-source  over 
single-source  information  described  in  [5], 

To  increase  the  realism  of  our  model,  we  introduce  the  notion  of  the  “ensembling” 
of  information  packages;  i.e.,  we  sequence  information  packages  for  processing  in 
ensembles — or  related  groups — rather  than  individually.  The  problem  of  interest,  then,  is 
to  develop  information  flow  models  for  the  case  where  information  packages  can  be 
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ensembled,  together  with  optimal  sequencing  results  and/or  heuristic  processing 
strategies.  We  introduce  two  new  concepts  here  that  represent  a  small  step  toward  that 
end. 

The  first  of  these  is  known  as  the  conditional  value  function.  As  we  observe 
above,  the  value  of  processing  a  particular  infonnation  package  is  dependent  upon  what 
has  already  been  processed.  We  define  a  conditional  value  function,  denoted  ValAj  (t) ,  as 

the  value  of  information  package  i  at  time  t  given  that  information  package  j  has  already 
been  processed.  It  would  be  safe  to  assume,  as  we  have  with  the  original  value  function, 
that  these  conditional  value  functions  are  piecewise  constant,  non-increasing  functions. 
Defining  these  functions  would  rely  heavily  on  existing  sensor  data  fusion  work  (see  [5], 
for  example). 

The  second,  related,  concept  is  that  Oval  2  can  be  thought  of  as  containing  subsets 
of  information  packages  rather  than  one  large  set.  Information  packages  within  subsets 
may  be  associated  by  the  targets  or  target  sets  they  describe,  the  geographical  region  they 
describe,  the  time  they  arrived,  or  other  factors.  As  such,  it  is  likely  that  subsets  within 
Oval  2  will  have  nonempty  intersections.  We  may  wish  to  consider  an  entire  subset  as  an 
ensemble  to  be  processed  together,  or  only  certain  portions  of  a  subset. 

Developing  an  infonnation  flow  model  that  treats  information  package  as 
ensembles  rather  than  individually  would  represent  an  important  achievement  in  that  it 
would  relate  two  previously  disjoint  modeling  areas:  dynamic  information  flow  and 
sensor  data  fusion. 

2.  Additional  Expected  Realized  Value  Results 

Let  Valt  (t)  =  (3,6  a,t  be  the  value  function  for  information  package 
i,  i  =  1,2 ,...,«.  Let  information  package  service  times  be  exponentially  distributed  with 
rate  //  and  define  the  permutation  mapping  <p(j)  to  be  the  information  package 
processed  z'th .  From  (3.40),  then,  we  know  that 
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(5.3) 


£[r]=iX)  n 

y=i  (=1 


This  expression  represents  the  expected  cumulative  realized  value  obtained  by  processing 
the  n  information  packages  in  any  sequence  cp .  The  following  questions  remain  open: 


•  Suppose  we  process  the  n  information  packages  in  the  sequence 
determined  by  heuristic  H5.  Does  a  closed  fonn  expression  for  E^Vr  J 
exist  and  can  it  be  derived? 

•  If  not,  can  we  define  reasonable  bounds  on  fs[F'J  for  n  information 
packages  processed  in  accordance  with  H5? 

•  Can  we  derive  exact  expressions  (or  bounds)  for  E  j^V  J  for  other  cases  of 

the  problem  (e.g.,  different  service  time  distributions  and  different  value 
function  types)? 


Answers  to  these  questions  would  provide  a  theoretical  basis  for  choosing  a 
heuristic  processing  strategy  to  employ  in  addition  to  the  empirical  evidence  shown  in 
Chapter  IV. 


3.  Additional  Theoretical  Results 

In  Chapter  III,  we  derived  results  for  simplified  cases  of  the  optimal  control 
problem  (2.4).  It  would  be  desirable  to  extend  these  results  to  at  least  the  following 
cases: 

•  Mixed  value  function  type  case  (suspected  to  be  NP-hard) 

•  Multi-processor  case  ( m  >2) 

4.  Extensions  to  Dynamic  Information  Flow  Modeling 

In  Chapter  II,  we  construct  a  dynamic  information  flow  model  that  describes  the 
flow  of  infonnation  between  Ovals  2  and  3  of  the  Dynamic  Model  of  Situated  Cognition. 
Figure  3  depicts  the  notion  of  developing  information  flow  models  for  other  portions  of 
the  DMSC  toward  the  ultimate  goal  of  tying  these  “sub-models”  together  to  fonn  a 
cohesive  representation  of  dynamic  information  flow  between  Ovals  1  and  6.  For 
example,  the  analysis  in  Chapter  III  is  all  based  on  the  assumption  that  the  number,  Q,  of 
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information  packages  to  be  sequenced  is  pre-determined,  or  fixed.  We  model  the  static 
scheduling  problem  as  if  the  commander  has  no  control  over  the  contents  of  Oval  2;  he 
only  controls  how  the  contents  of  Oval  2  are  scheduled  for  processing.  This  is  not  the 
case  in  actuality.  The  commander  controls  his  organic  sensor  assets;  he  determines  how 
they  are  allocated,  which  directly  affects  the  number  of  information  packages  present  in 
Oval  2  at  any  time,  and  even  the  type  of  IPs  present.  So,  a  model  is  required  that  relates 
the  contents  of  Oval  2  to  a  commander’s  sensor  allocation  strategy  (the  control).  This 
model  could  be  thought  of  as  “residing”  between  Ovals  1  and  2  of  the  dynamic  model  of 
situated  cognition,  describing  how  information  flows  between  “ground  truth”  (Oval  1) 
and  the  set  of  information  collected  (Oval  2).  Tying  such  a  model  together  with  the 
dynamic  information  flow  model  presented  in  Chapters  II  and  III  would  represent  a 
significant  step  toward  the  ultimate  goal  of  mathematically  modeling  the  flow  of 
information  between  Ovals  1  and  6. 

5.  Implementing  a  “Warm”  Start  in  the  Dynamic  Scheduling  Simulation 

The  dynamic  scheduling  problem  we  have  investigated  assumes  a  “cold”  start; 
i.e.,  zero  information  packages  are  present  in  the  system  at  time  zero.  Another  realistic 
starting  condition  would  be  a  “warm”  start  where  the  scheduling  process  does  not  begin 
until  some  time  after  information  packages  have  begun  arriving.  Preliminary  findings 
seem  to  indicate  that  there  is  no  significant  difference  between  the  two  approaches,  but  a 
more  comprehensive  simulation  investigation  is  required  before  drawing  any  firm 
conclusions. 
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APPENDIX  A.  SIMULATION  CODE  (MATLAB)  AND 
ADDITIONAL  EXPERIMENTAL  RESULTS 


A.  STATIC  SIMULATION  CODE 

%  StaticPolicyTest .m 

%  compares  svc  policies  for  sets  of  IPs 

%  assumes  a  fixed  #  of  IPs  to  be  sequenced  with  no  new  arrivals 
%  Created  10  Feb  2008. 


%  Assumptions: 

%  -  exponential,  linear,  or  step  value  functions  (arbitrary  decay  rate 
and  initial  value) 

o, 

o 

%  No  new  arrivals 

%  3  "classes"  of  IPs  -  determines  service  time  distributions 
clear  all; 

%  define  parameter  matrix  -  each  row  represents  a  run. 

%  cols  are:  run  #,  reps,  Q,  ValFctMix,  k  alpha,  theta  alpha,  k  beta, 
theta_beta, 

%  k_tau,  theta_tau,  k^step,  theta^step,  E_S1,  E_S2,  E_S3.  Tau  is  the 
ending  time  for  step 

%  functions.  k  step  and  theta_step  are  the  parameters  of  the  gamma 
%  distribution  from  which  the  initial  value  for  step  value  functions  is 
%  drawn. 

%  ValFctMix  is  1  if  all  exp  val  fcts;  2  if  all  linear  val  fcts; 

%  3  if  all  step  fcts;  4  if  "mix"  -  equal  mix  for  now. 
param=importdata ( ' runmatrix . csv ' ) ; 

overall=zeros (size (param, 1) , 6) ; 
tic 

for  run=l : size (param, 1 ) 

reps=param (run, 2 ) ;  %#  of  iterations 
Q=param (run, 3) ;  %  number  of  IPs; 

ValFctMix  =  param ( run, 4 ) ; 

results=zeros (reps, 6) ;  %  tracks  output  for  each  run 

count= [0  0  0 ] ; 
for  z=l : reps 

count= [count ( 1 ) +1  run  size (param, 1 ) ] 

%rand (' twister ',  5489)  ;  %  returns  the  same  values  each  time 
rand (' twister ' ,  sum ( 1 00*clock) ) ;  %  random  seed  driven  by  clock 
time;  returns  different  values  each  time 
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User  defined  inputs: 


numclass  =  3; 

E  S= [param (run, 13) ;  param (run, 14 ) ;  param (run, 15) ] ;  %  mu  i  are 
mean  service  times  for  IPs  of  class  i  (service  time  distributions  are 
assumed  to  be  exponential) 

alpha  =  gamrnd (param (run, 5) , param (run, 6) , Q, 1) ;  %  Q-long  vector 
of  decay  rates  (lin  and  exp) 

beta  =  gamrnd (param ( run, 7 ), param ( run, 8 ), Q, 1 ) ;  %  Q-long  vector 
of  initial  values  (lin,  exp,  and  step) 

%step  val  =  gamrnd (param ( run, 1 1 ), param (run, 12 ), Q, 1 ) ;  %  vector 
of  constant  vals  for  step  val  fcts 

tau  =  gamrnd (param ( run, 9 ), param ( run, 1 0 ), Q, 1 ) ;  %  vector  of 
ending  times  (step  val  fcts  only) 

class  =  ceil (numclass*rand (Q, 1 ) ) ;  %  Q-long  vector  of  IP  classes 

expectedsvctime  =  E_S (class); 

ValFctForm=zeros (Q, 1) ; 
if  ValFctMix  ==  1 

ValFctForm  =  ones(Q,l);  %if  all  exp  val  fcts 
elseif  ValFctMix  ==  2 

ValFctForm  =  2*ones(Q,l);  %  if  all  linear  val  fcts 
elseif  ValFctMix  ==  3 

ValFctForm  =  3*ones(Q,l);  %  if  all  step  val  fcts 

else 

%ValFctForm  =  ceil ( 3*rand (Q, 1 ) ) ;  %  creates  a  uniform  mix  of 
all  three  val  fct  types. 

vr=rand (Q, 1 ) ; 

ValFctForm  =  (vr< . 75) +2* (vr>= . 75  &  vr< . 9 ) +3* (vr>= . 9) ; 

end 


%initialize  the  IP  characteristic  data  array 
IP=zeros (Q, 12 ) ; 

IP ( : , 1)  =  1:Q; 

IP ( : , 2 )  =  alpha; 

IP ( : , 3)  =  beta; 

IP ( : , 4 )  =  class; 

IP(:,5)  =  expectedsvctime;  %  this  is  E[S_sub_q] 
%IP(:,6)  =  step_val; 

IP ( : , 7 )  =  tau; 

IP ( : , 8 )  =  ValFctForm; 


%  col's  9-12  of  IP  are  the  vector  "param"  for  input  to  the 
%  function  ValueFunction . 

for  j  j  =1 : Q 

if  IP(jj,8)==l  %  exp  val  fct 
IP (j j , 9) =IP ( j j , 2) ;  %  alpha 
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end 


IP ( j j , 10) =IP ( j j , 3) ;  %  beta 
elseif  IP(jj,8)==2  %  lin  val  fct 
IP ( j j , 9) =IP ( j j , 2) ;  %  alpha 
IP ( j j , 10) =IP ( j j , 3) ;  %  beta 
elseif  IP(jj,8)==3  %  step  val  fct 
IP(jj,9)=IP(jj,3)  ; 

IP ( j j , 10) =IP  ( j j , 7 ) ; 

end 


param  val=IP ( : , 9 : 12) ; 
tic 

if  Q<=10  %only  find  exact  solution  for  small  Q 

seq  =  single (perms ( 1 : Q) ) ;  %a  Q!  x  Q  matrix  of  all  possible 
permutations  of  1:Q;  each  row  represents  a  processing  sequence. 

seq  =  single ([seq  zeros ( factorial (Q) , 1 )]) ;  %  appends  a 
column  of  zeros  at  the  end  of  seq  (overwrites  the  original  to  save 
memory) 

%  this  last  column  will  hold  realized 
%  values  for  each  row  (sequence  of  IPs) 

max  val  =0;  %  used  to  track  optimal  value  -  updated  while 
simulation  runs 

for  j =1 : factorial (Q) 

tl=IP (seq ( j , 1 :Q) , 5) ; 

t  =  cumsum(tl);  %  vector  of  expected  server  availble 

times 

form=IP (seq ( j , 1 :Q) , 8) ;  %  val  fct  forms 
a=zeros (Q, 1 ) ;  %arrival  time  is  zero 
param_shuf f le=param_val ( seq ( j  ,  1 :  Q)  ,  : )  ; 

val  k  =  diag (ValueFunction (t,  form,  a,  param  shuffle)); 
realized  val  =  sumfval  k) ; 

%  for  k=l : Q 

%  t  =  t  +  IP (seq ( j , k) , 5) ;  %  updates  running  expected 

time 

%  val_k  =  ValueFunction (t,  IP (seq ( j , k) , 8) ,  0, 
param_val ( seq ( j  ,  k)  ,  :  )  )  ; 

%  realized  val  =  realized  val  +  val  k; 

%end 

seq(j,Q+l)  =  realized_val; 

if  realized  val  >  max  val 

max  val  =  realized  val;  %  keeps  track  of  current 

optimal  value 

end 


end 

OptimalValue  =  max  val; 
results (z,l)  =  OptimalValue; 
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end 

toe 


%HEURISTIC  #1  DOES  NOT  APPLY  TO  STEP  FUNCTIONS  -  SO  WE  OMIT. 

%  Heuristic  #2:  Really  Greedy  (sorts  on  value  at  t  =  0, 
without  resorting) 

%  so,  really,  this  just  sorts  on  Beta. 

IPh2  =  sortrows (IP, -3) ;  %  sorts  by  Beta  in  descending  order 
param  val=IPh2 ( : , 9 : 12 ) ; 

t=0; 
rvh2=0 ; 
for  k=l : Q 

t  =  t  +  IPh2(k,5);  %  updates  running  expected  time 
val  k  =  IPh2(k,3)  *  exp (-IPh2 (k, 2) *t) ; 

val  k  =  ValueFunction (t,  IPh2(k,  8),  0,  param  val(k,  :)); 
rvh2  =  rvh2  +  val  k; 

end 

results  (z, 3)  =  rvh2; 

%  Heuristic  #3:  Pretty  Greedy  (same  as  #2,  only  this  time  we 
re-sort  after 

%  every  IP  selection,  evaluating  the  value  function  at  current 

time) . 

%  Expect  performance  to  be  no  worse  than  h2 . 

[r  c]  =  size (IP) ; 

IPh3=zeros (r, c+1 ) ; 

IPtemp  =  [IP  zeros (r,l)]; 

t=0; 

for  k=l : Q 

param  val=IPtemp ( : , 9 : 12 ) ; 
a= zeros (size (IPtemp, 1) , 1) ; 

IPtemp (:, c+1 )= (ValueFunction (t,  IPtemp (:, 8),  a, 
param  val)) ' ;  ^current  value 

IPtemp  =  sortrows ( IPtemp, - (c+1 )) ; 

IPh3(k, :)  =  IPtemp(l,:); 

t  =  t  +  IPtemp (1,5);  %  updates  the  clock 

IPtemp (1,:)  =  [];  %deletes  first  row  of  IPtemp  (the  row 
that  was  just  copied  to  IPh3) 

%  note  -  IPtemp  loses  a  row  each  iteration  - 
%  it  eventually  disappears! 


end 

t=0; 
rvh3=0 ; 

param  val=IPh3  ( : , 9  : 12 )  ; 
for  k=l : Q 
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t  =  t  +  IPh3(k,5);  %  updates  running  expected  time 
val  k  =  ValueFunction (t,  IPh3(k,  8),  0,  param  val(k,  :)); 
rvh3  =  rvh3  +  val  k; 

end 

results  (z, 4)  =  rvh3; 

%  Heuristic  #4:  Somewhat  Greedy  -  look  ahead  to  predict  value 

AFTER 

%  service  time;  resort  each  time 

IPh4=zeros (r,c+l) ; 

IPtemp  =  [IP  zeros (r,l)]; 
t=0; 

for  k=l:Q 

param  val=IPtemp ( : , 9 : 12) ; 
a=zeros (size (IPtemp, 1) , 1)  ; 

IPtemp ( : , c+1) =diag (ValueFunction (t+IPtemp ( : , 5) , 

IPtemp (:, 8),  a,  param_val) ) ;  %  val  at  end  of  svc 
IPtemp  =  sortrows (IPtemp, - (c+1) ) ; 

IPh4(k, :)  =  IPtemp(l,:); 

t  =  t  +  IPtemp (1,5);  %  updates  the  clock 

IPtemp (1, :)  =  [];  %deletes  first  row  of  IPtemp  (the  row 
that  was  just  copied  to  IPh3) 

%  note  -  IPtemp  loses  a  row  each  iteration  - 
%  it  eventually  disappears! 

end 

t=0; 
rvh4=0 ; 

param  val=IPh4  ( : , 9  : 12 )  ; 
for  k=l:Q 

t  =  t  +  IPh4(k,5);  %  updates  running  expected  time 
val  k  =  ValueFunction (t,  IPh4(k, 8),  0,  param  val(k, :)); 

rvh4  =  rvh4  +  val  k; 

end 

results  (z, 5)  =  rvh4; 

%Heuristic  #5:  "bang  for  buck"  -  descending  sort  on 
%Val_q (E [d_q] ) /E [ S_q]  -  resort  each  time 

IPh5=zeros (r,c+l) ; 

IPtemp  =  [IP  zeros (r,l)]; 

t=0; 

for  k=l : Q 

param  val=IPtemp ( : , 9 : 12 ) ; 
a= zeros (size (IPtemp, 1) , 1) ; 

IPtemp ( : , c+1) = (diag (ValueFunction (t+IPtemp ( : , 5) , 

IPtemp (:, 8),  a,  param_val) )). /IPtemp (:, 5) ;  %  val  at  end  of  svc 
IPtemp  =  sortrows(IPtemp,-(c+l)); 

IPh5(k, :)  =  IPtemp(l,:); 
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t  =  t  +  IPtemp(l,5);  %  updates  the  clock 
IPtemp(l, :)  =  [];  %deletes  first  row  of  IPtemp  (the  row 
that  was  just  copied  to  IPh3) 

%  note  -  IPtemp  loses  a  row  each  iteration  - 
%  it  eventually  disappears! 

end 

t=0; 
rvh5=0 ; 

param  val=IPh5  ( : , 9  : 12)  ; 
for  k=l : Q 

t  =  t  +  IPh5(k,5);  %  updates  running  expected  time 
val  k  =  ValueFunction (t,  IPh5(k, 8),  0,  param  val(k, :)); 
rvh5  =  rvh5  +  val  k; 

end 

results (z, 6)  =  rvh5; 


end 

%Normalize  results  matrix  -  either  to  exact  optimal  (if  Q<=10)  or 
to  rvh3 

%(the  most  reasonable  greedy  heuristic),  then  construct  numerical 
%histogram. 

[row, col]  =  size (results) ; 
results  norm=zeros (row,  col)  ; 

if  Q<=10  %  we  normalize  on  exact  optimal  results  in  this  case 
for  i=l : col 

results_norm (:, i ) =results (:, i ). /results (:, 1 ) ;  %  first 
column  of  results  is  optimal 
end 

bin= (0 : . 05 : 1) ' ;  %bin  sizes  for  histogram 
h  =  hist (results  norm, bin) ; 
overall (run, 1 ) =param (run, 1 ) ; 
for  bb=2 : 6 

overall (run, bb) =sum (results ( : , bb) ) /sum (results  (:, 1 )) ; 

end 

end 

if  Q>10  %  we  normalize  on  rvh3  in  this  case 
for  i=l : col 

results_norm (:, i ) =results (:, i ). /results (:, 4 ) ;  %  fourth 
column  of  results  is  rvh3 
end 

bin= ( 0 : .5:100) ' ; 
h  =  hist (results  norm, bin); 


end 


h= [bin  h] ; 

csvwrite ( [ ' hist  run  ' , num2str (run) , ' . csv ' ] , h) ; 
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csvwrite ( [ ' norm  results  run  num2str (param (run, 1 csv '], results  nor 

m)  ; 


%bar3 (bin,  h) 


end 

csvwrite ( ' overall . csv ' , overall ) ; 
toe 


B.  VALUE  FUNCTION  CODE 

This  function  is  incorporated  as  a  sub-routine  in  both  the  static  and  dynamic 
simulation  codes. 


%%  y  =  ValueFunction (t,  ValFctForm, a, param) 

%%  This  m-file  outputs  Val  j (t)  for  independent  variable  t  (scalar  or 
%%  vector) .  ValFctForm  indicates  the  type  of  value  function  (1  = 

%%  exponential;  2  =  linear;  3  =  step) .  a  is  the  arrival  time  for 
package 

%%  "j"  -  the  given  package.  param  is  a  row  vector  containing  function 
%%  parameters . 

%%  assume  t  is  k  x  1;  ValFctForm  and  a  are  n  x  1;  param  is  n  x  4 

%%  Version  2,  27  June  2007 
%%  Donovan  Phillips 

function  y  =  ValueFunction (t,  ValFctForm, a, param) 

%  convert  all  input  to  column  vectors  if  not  already 

if  (size (t, 1) ==1)  &&  (size (t, 2) ~=1) 

t=t '  ; 

end 

if  (size (ValFctForm, 1) ==1)  &&  (size (ValFctForm, 2) ~=1) 

ValFctForm=ValFctForm ' ; 

end 

if  (size (a, 1) ==1)  &&  (size (a, 2) ~=1) 

a=a '  ; 

end 

if  size (param, 1 )  ~=  max (size (a)) 

param=param' ; 

end 
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vl  =  ValFctForm  ~=1;  v2  =  ValFctForm  ~=  2;  v3  =  ValFctForm  ~=  3; 

v  =  vl  . *  v2  .*  v3;  %  v  will  contain  non-zero  entries  if  ValFctForm 
contains  entries  other  than  1,  2,  or  3 . 

if  sum(v)  >  0 

error ( [num2str (ValFctForm) , '  is  an  invalid  functional  form.  Must 
use  1 ,  2 ,  or  3 . ' ] ) 
end 

if  (size (ValFctForm, 1)  ~=  size (a, 1))  | |  (size (ValFctForm, 1)  ~= 

size (param, 1) ) 

error (' the  vectors  ValFctForm,  a,  and  param  must  have  the  same 
number  of  rows ' ) 
end 


size  param  =  size (param, 2 ) ; 
if  size  param  ~=  4 

error ('param  must  have  4  elements  for  every  package') 

end 

size_t  =  max (size (t) ) ;  size_a  =  size (a, 1); 

%  organize  input  data 


j  =  cumsum (ones ( size_a, 1 ) ) ; 
input  =  [j  ValFctForm  a  param] ; 

input  sort  =  sortrows (input, 2) ;  %  sorts  on  ValFctForm  in  ascending 
order 

il  =  input_sort ( : , 2 )  ==  1;  i2  =  input_sort ( : , 2 )  ==  2;  i3  = 
input_sort ( : , 2 )  ==  3; 


[cl  hi]  =  min(il);  %  cl  is  the  min  value,  hi  is  the  index  of  the  first 
instance  of  this  min  value 

%  so,  row  hi  of  input  sort  is  the  first  row  having 
%  ValFctForm  other  than  1  (note:  could  be  the 

first 


%  row  -  ie,  there  ARE  no  1 
[c2  h2]  =  max(i2);  %  if  c2=l,  then  there  is  at 
the  first  of  these  has  index  h2 
[c3  h3]  =  max(i3);  %  if  c3  =  1  (usually  will), 
first 

%  "3"  in  ValFctForm. 

sumil=sum ( il ) ;  sumi2=sum ( i2 ) ;  sumi3=sum ( i3 ) ;  % 
1 ' s  in  each  vector 


s) 

least  one  ValFctForm=2 ; 
h3  is  the  index  of  the 

indicates  the  number  of 


%  Note:  if  hi  =  h3,  then  there  are  no  "2's"  in  ValFctForm;  if  hl=l, 
there 

%  are  no  l's;  -  need  to  cover  all  these  bases... 

%  initialize  output  matrix,  y 
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y  =  zeros (size  t,  size  a);  %  output  will  have  a  column  for  every  pkg 
input 


-k  ~k  -k  ~k  -k  -k  -k  ★  -k  -k  -k  -k  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  9c  -k  ~k  -k  -k  -k  -k  -k  -k  -k  ~k 

%  Exponential  Value  Function  Form  (ValFctForm  =  1) 

-k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  ~k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  -k  ~k  -k  -k 


if  sumil  >  0  %  ie,  if  there  is  at  least  one  row  of  input  sort  with 
ValFctForm=l 

input  sortl  =  input  sort ( 1 : sumil ;  %  only  the  rows  with 
ValFctForm=l 

for  k=l : sumil 

pkg  =  input_sortl ( k, 1 ) ; 

a  =  input  sortl (k,  3);  %arrival  time  of  the  kth  pkg 
pi  =  input  sortl (k, 4);  p2  =  input  sortl (k, 5);  %  function 
parameters 

arriv  =  t  >=  a;  %  arriv  is  a  vector  -  same  length  as  t  - 

%  that  is  0  if  t<arrival  time 


1 


out  =  p2*exp (-pi* (t-a) 
nonnegative) 

y(:,pkg)  =  arriv. *out; 


%  otherwise 

;  %  pi  =  alpha,  p2  =  beta  (both 
%  only  considers  values  of  t  greater 


than  arrival  time 
end 


end 


SZ-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Linear  Value  Function  Form  (ValFctForm  =  2) 

SZ-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


if  sumi2  >  0  %  ie,  there  is  at  least  one  ValFctForm  =  2 

input  sort2  =  input  sort (h2 : (h2+sumi2-l ) , : ) ;  %  only  the  rows  of 
input  sort  having  ValFctForm=2 
for  k  =  1 : sumi2 

pkg  =  input_sort2 ( k, 1 ) ; 

a  =  input  sort2(k,3);  %arrival  time  of  the  kth  pkg 
pi  =  input  sort2(k,4);  p2  =  input  sort2(k,5);  %  function 
parameters 

arriv  =  t  >=  a;  %  arriv  is  a  vector  -  same  length  as  t  - 

%  that  is  0  if  t<arrival  time 


1 


%  otherwise 

out  =  max (0,  p2  -  pi* (t-a) ) ;  %  pl=alpha,  p2=beta, 
nonnegative 

y(:,pkg)  =  arriv. *out;  %  only  considers  values  of 
than  arrival  time 
end 


both 

t  greater 


end 


Q~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  Step  Value  Function  Form  (ValFctForm  =  3) 


97 


S-'k  k  ~k  ~k  ~k  ~k  ~k  k  ~k  k  ~k  k  ~k  ~k  ~k  ~k  ~k  k  ~k  ~k  ~k  ~k  ~k  ~k  ~k  k  -k  k  ~k  k  ~k  ~k  ~k  -k  ~k  ~k  ~k  ~k  ~k  k  ~k  k  ~k  ~k  ~k  -k  ~k  ~k  ~k 


if  sumi3  >  0  %  ie,  there  is  at  least  one  ValFctForm  =  3 

input  sort3  =  input  sort (h3 : (h3+sumi3-l ) , : ) ;  %  only  the  rows  of 
input  sort  having  ValFctForm=3 
for  k  =  1 : sumi3 

pkg  =  input_sort3 ( k, 1 ) ; 

a  =  input  sort3(k,3);  %arrival  time  of  the  kth  pkg 
pi  =  input_sort3 ( k, 4 ) ; 
p2  =  input_sort3 ( k, 5 ) ; 
p3  =  input_sort3 (k, 6) ; 

p4  =  input  sort3(k,7);  %  function  parameters 

%arriv  =  t  >=  a;  %  arriv  is  a  vector  -  same  length  as  t  - 

%  that  is  0  if  t<arrival  time 


1 


%  otherwise 

a_to_p2  =  (t>=a)  &  (t<p2);  %  true  if  a<=t<p2 

p2_to_p4  =  (t>=p2)  &  (t<p4); 


y(:,pkg)  =  pl*a_to_p2  +  p3*p2_to_p4; 

end 

end 


C.  DYNAMIC  SIMULATION  CODE 

%  simulation  to  test  various  service  (processing)  policies 

%  first,  must  create  information  packages,  with  associated  arrival 
times 

%  (initially  assumed  exponential),  service  times  (again,  exponential 
%  different  rates  for  each  class  of  info  pkg),  and  value  function. 

clear  all; 
tic 

parameters=importdata ( ' Dynamiclnput . csv ' ) ; 

NumRuns  =  size (parameters, 1) ; 

overall  mean=zeros (NumRuns, 8) ;  %  first  col  is  run  #;  then  mean/stdev 
for  h5,  h4,  h3,  LIFO,  Random,  StatusQuo,  &  FIFO 
overall  std=zeros (NumRuns , 8 ) ; 

for  run  =  1: NumRuns 

NumReps  =  parameters (run,  2) ; 

RunResults=zeros (NumReps+2 , 7 ) ;  %  last  two  rows  hold  mean,  stdev; 
cols  are  h5  -  FIFO  (same  order  as  above) 
for  rep  =  1: NumReps 

run 

rep  %visual  status  output 

%rand (' twister ', 5489) ;  %  returns  the  same  values  each  time 
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rand (' twister ' ,  sum ( 1 00*clock) ) ;  %  random  seed  driven  by  clock 
time;  returns  different  values  each  time 

%  USER  INPUT  SECTION  ****************** 

NumClass  =  parameters (run, 4 ) ;  %  k  =  1...3  classes  (1  - 

"EO/IR";  2  -  "Still  Imagery";  3  -  "Voice") 

NumValFctForm  =  parameters (run, 27) ;  %  j  =  1...3  ValFctForm  (1 

=  exp;  2  =  linear;  3  =  step) 

ValFctMix  =  parameters (run, 6 : 8) ;  %  percent  of  val  fcts  that  are 
exp,  lin,  and  step,  respectively. 

%  Number  of  elements  must  equal  NumValFctForm. 

%  elements  must  sum  to  1 . 

NumPkg  =  parameters (run, 3) ;  %  j  =  1...1000  pkgs 
NumSvr  =  parameters (run, 5) ;  %  number  of  information  package 
processors 

lambda  =  parameters (run, 21 : 23) ';  %  arrival  rates  for  each  class 
mu  =  parameters (run, 24 : 26) ';  %  service  RATES  by  class 
gamparam  = 

[parameters ( run, 9:10) ; parameters ( run, 11:12) ; parameters ( run, 13:14) ; par am 
eter s ( run ,15:16); parameters ( run ,17:18); parameters ( run ,19:20)]; 

%  gamparam  is  a  2-column  matrix  of  parameters  for  gamma 
distributions  from 

%  which  the  parameters  for  each  value  function  type  are  drawn. 

The  two 

%  cols  represent  "k"  and  "theta"  -  the  two  parameters  of  the 
gamma  distn; 

%  each  row  represents  a  parameter  as  follows:  1st  and  2nd  rows 
represent 

%  alpha  and  beta,  respectively,  for  exponential  value 
functions;  3rd  and 

%  4th  rows  -  alpha  and  beta  for  linear  val  fcts;  5th  row  is 
beta  (initial 

%  value)  for  step  fcts;  6th  row  is  "tau"  -  ending  time  for  step 
val  fcts. 


%  val  fcts  are  defined  as  follows.  Exp:  val 


%  END  USER  INPUT  SECTION  ***************** 

%  generate  class  and  arrival  times  (each  of  length  NumPkg) : 

mean  =  1. /lambda; 

inter  arrl  =  exprnd (mean ( 1 ), NumPkg,  1 ) ;  arrl  = 
cumsum( inter  arrl); 

inter  arr2  =  exprnd (mean (2 ), NumPkg, 1 ) ;  arr2  = 
cumsum( inter  arr2); 

inter  arr3  =  exprnd (mean ( 3 ), NumPkg, 1 ) ;  arr3  = 
cumsum( inter  arr3); 

arrlc  =  [arrl  ones (NumPkg, 1 )] ; 

arr2c  =  [arr2  2*ones (NumPkg, 1 )] ; 

arr3c  =  [arr3  3*ones (NumPkg, 1 )] ; 
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arrc  =  [arrlc;  arr2c;  arr3c] ; 
arr  sort  =  sortrows (arrc, 1 ) ; 


3) 

last 


class  =  arr  sort ( 1 : NumPkg, 2 ) ;  %  information  pkg  class  (1,  2  or 

a  =  arr  sort ( 1 : NumPkg, 1 ) ;  %  arrival  times  -  sorted  first  to 

ExpSvcTime  =  l./mu;  %  expected  svc  times  for  each  class 
es  =  ExpSvcTime (class) ;  %  expected  svc  times  for  each  IP 
endtime  =  a (NumPkg) *1.01; 


%  generate  ValFctForm  vector  according  to  ValFctMix: 


vr=rand (NumPkg, 1 ) ; 

ValFctForm=zeros (NumPkg,  1) ; 
vt  =  [0  cumsum (ValFctMix) ] ; 
for  j =1 : NumValFctForm 

ValFctForm=ValFctForm  +  j*(vr  >  vt(j)  &  vr  <=  vt(j+l)); 

end 


%  generate  service  times 
s  =  exprnd (1 . /mu (class) ) ; 

%  initialize  and  construct  the  parameter  array: 

param  =  zeros (NumPkg, 4 ) ;  %keep  the  4  cols:  ValueFunction . m 
requires  it  (cols  3  &  4  will  be  zeros) 
for  i=l: NumPkg 

if  ValFctForm ( i ) ==1  %  exponential 

param(i,l)  =  gamrnd (gamparam ( 1 , 1 ) , gamparam ( 1 , 2 ) ) ;  % 

alpha 

par am (i, 2)  =  gamrnd (gamparam (2 , 1 ), gamparam (2 , 2 )) ;  % 

beta 

elseif  ValFctForm ( i )  ==  2  %  linear 

param(i,l)  =  gamrnd (gamparam ( 3 , 1 ), gamparam ( 3 , 2 )) ;  % 

alpha 

par am (i, 2)  =  gamrnd (gamparam ( 4 , 1 ), gamparam ( 4 , 2 )) ;  % 

beta 

elseif  ValFctForm ( i )  ==  3  %  step 

param(i,l)  =  gamrnd (gamparam ( 5 , 1 ), gamparam ( 5 , 2 )) ;  % 

beta 

param(i,2)  =  a(i)  + 

gamrnd (gamparam ( 6, 1 ), gamparam ( 6, 2 )) ;  %  tau 
end 

end 


%  determine  "realized  value"  over  time  for  various 
%  service  policies: 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 
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%  SERVICE  POLICY  #5:  (value  at  expected  svc 

completion) / (expected  svc 
%  time)  H5 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

IPs5  =  [ ( 1 : NumPkg) '  a  zeros (NumPkg, 1 )] ;  %  first  col  is  IP 
permanent  index 

%  sequenced  by  arrival  time;  second  col  is  arrival  time;  third 
col  is  an 

%  indicator  which  is  1  if  the  IP  has  already  been  sent  to 
server  or  0  if 
%  not . 

SvrWkld5  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time 
assigned  to  each  server. 

SvrAsgn5  =  zeros (NumPkg, 3 ) ;  %each  row  represents  a  package 
(sorted  by  arrival 

%  time,  just  like  IPs);  first  column 
%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time. 
ct5  =  zeros (NumPkg, 1 ) ; 

ct5(l)  =  a(l);  %  ct="clock  time"  -  starts  upon  first  arrival. 

%  note  -  ct  is  the  time  an  assignment  decision  is  made  -  at 
%  least  one  server  must  be  available. 

for  i  =  1: NumPkg 

%  extract  from  IPs  those  pkgs  that  have  1)  arrived  (i.e., 
a<=ct)  and  2) 

%  those  that  have  not  yet  been  sent  to  svr  (3rd  col  of  IPs 

is  zero) 


IPtemp  =  sortrows (IPs5, 2) ;  %sorts  in  ascending  order  on 

arrival  time 

arr  =  IPtemp (:, 2)  <=  ct5(i);  %  arr  =  1  if  IP  has  arrived 
(a<=ct5) ;  arr  =  0  otherwise 

%  note:  arr  is  always  size  (NumPkg  x  1) . 

%  Also,  sum(arr)  =  #  of  arrivals  as  of 
%  ct  ( i )  . 

[cc  dd]  =  min (arr);  %  dd  is  the  index  of  the  first  zero 
term,  if  there  is  a  zero 

if  cc  ~=1  %  ie,  if  there  is  at  least  one  IP  that  hasn't 

"arrived" 

rr=size (IPtemp, 1) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp (dd: rr, : )  =  [];  %  deletes  all  rows  w/arrival 
times  greater  than  current  time 
end 

IPtemp  =  sortrows ( IPtemp,  -3);  %  sorts  in  descending  order 
on  col  3  (sent  to  svr) 

[aa  bb]  =  min (IPtemp (:, 3) ) ;  %  bb  is  the  row  containing  the 

first  zero 
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IPtemp (1 : (bb-1) , : )  =  [];  %  deletes  rows  having  svr  assg 
indicator  of  1  (sent  to  svr) 


current 
IPtemp 
each  IP 


r  =  size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 
%  here ' s  where  we  implement  h5 : 

IPtemp  =  [IPtemp  zeros (r,  1) ] ;  %  add  a  column  to  hold  the 
value 

rows  =  IPtemp (:,1);  %  these  are  the  IP  indices  remaining  in 

ExpDep  =  ct5(i)  +  es(rows);  %  expected  departure  time  for 
if  it  is  sent  to  svc  at  time=ct(i) 


IPtemp ( : , 4 )  = 

(diag (ValueFunction (ExpDep, ValFct Form ( rows ) , a (rows ) , par am ( rows ,:)))). /e 
s (rows) ; 

%  computes  h5  for  each  IP  left  in  IPtemp 

IPtemp  =  sortrows ( IPtemp, -4 ) ;  %  sorts  on  h5,  descending 

IPpick  =  IPtemp (1,1);  %choose  IP  w/highest  h5  value 

IPs5 (IPpick, 3)  =  1;  %IP  w/highest  h5  value  goes  to  server 

[v  index]  =  min ( SvrWkld5 ) ;  %  index  is  the  next  avail  svr;  v 
is  the  time  it's  avail. 

SvrAsgn5 ( IPpick, 1 )  =  index;  %assigns  IPpick  to  the  next 
available  server 

if  a (IPpick)  <=  v 

SvrAsgn5 ( IPpick, 2 )  =  v; 

else 

SvrAsgn5 ( IPpick, 2 )  =  a (IPpick);  %  svc  commencement  is 
the  later  of  arrival  time  and  svr  avail  time. 

end 

SvrAsgn5 (IPpick, 3)  =  SvrAsgn5 ( IPpick, 2 )  +  s (IPpick);  %  svc 

completion  time 

SvrWkld5 ( index)  =  SvrAsgn5 ( IPpick, 3 ) ;  %  time  when  this 
server  is  next  available. 


%  update  clock  time  (except  for  the  last  iteration) : 
if  i  <  NumPkg 


if  (sum(arr) 
ct5  ( i  +  1 ) 

else 

ct5 ( i+1 ) 

end 


>=  i  +  1 )  &&  (min ( SvrWkld5 )  <  a (i+1)) 

=  min ( SvrWkld5 ) ; 

=  max (a (i+1),  min ( SvrWkld5 ) ) ; 


end 


end 

%  save  this  for  output  and  analysis  as  necessary. 
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col  1  is  IP 


IPdata5  =  [IPs5(:,l:2)  SvrAsgn5  zeros (NumPkg, 1 )] ;  % 
index  (sequenced  by  arrival  time) ; 

%  col  2  is  arrival  time;  col  3  is  server 
%  assigned;  col  4  svc  start  time;  col  5  is  svc 
%  completion  time;  col  6  is  realized  value. 

IPdata5(:,6)  =  diag (ValueFunction (IPdata5 ( : , 5) ,  ValFctForm,  a, 

param) ) ; 


IPdata5sort  =  sortrows (IPdata5, 5) ;  %  sort  on  svc  completion 

time 

IPdata5sort  =  [IPdata5sort  cumsum ( IPdata5sort ( : , 6 ) ) ] ;  %  add 
additional  col  for  cumulative  realized  value 

cutoff  =  IPdata5sort ( : , 5)  <=  endtime; 

[g2  h2 ] =min (cutoff ) ; 

RunResults (rep, 1 )  =  IPdata5sort (h2 , 7 ) ;  %  this  is  total  realized 
value  for  h5  for  this  rep 

^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  SERVICE  POLICY  #4:  Value  at  expected  svc  completion  h4 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

IPs4  =  [( 1 : NumPkg) '  a  zeros (NumPkg, 1 )] ;  %  first  col  is  IP 
permanent  index 

%  sequenced  by  arrival  time;  second  col  is  arrival  time;  third 
col  is  an 

%  indicator  which  is  1  if  the  IP  has  already  been  sent  to 
server  or  0  if 
%  not . 

SvrWkld4  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time 
assigned  to  each  server. 

SvrAsgn4  =  zeros (NumPkg, 3 ) ;  %each  row  represents  a  package 
(sorted  by  arrival 

%  time,  just  like  IPs);  first  column 
%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time. 
ct4  =  zeros (NumPkg, 1 ) ; 

ct4(l)  =  a(l);  %  ct="clock  time"  -  starts  upon  first  arrival. 

%  note  -  ct  is  the  time  an  assignment  decision  is  made  -  at 
%  least  one  server  must  be  available. 

for  i  =  1: NumPkg 

%  extract  from  IPs  those  pkgs  that  have  1)  arrived  (i.e., 
a<=ct)  and  2) 

%  those  that  have  not  yet  been  sent  to  svr  (3rd  col  of  IPs 

is  zero) 


IPtemp  =  sortrows ( IPs4 , 2 ) ;  %sorts  in  ascending  order  on 

arrival  time 
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arr  =  IPtemp(:,2)  <=  ct4(i);  %  arr  =  1  if  IP  has  arrived 
(a<=ct4);  arr  =  0  otherwise 

%  note:  arr  is  always  size  (NumPkg  x  1) . 

%  Also,  sum(arr)  =  #  of  arrivals  as  of 
%  ct ( i ) . 

[cc  dd]  =  min (arr);  %  dd  is  the  index  of  the  first  zero 
term,  if  there  is  a  zero 

if  cc  ~=1  %  ie,  if  there  is  at  least  one  IP  that  hasn't 

"arrived" 

rr=size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp (dd: rr, : )  =  [];  %  deletes  all  rows  w/arrival 
times  greater  than  current  time 
end 

IPtemp  =  sortrows ( IPtemp,  -3);  %  sorts  in  descending  order 
on  col  3  (sent  to  svr) 

[aa  bb]  =  min (IPtemp (:, 3) ) ;  %  bb  is  the  row  containing  the 

first  zero 


IPtemp ( 1 : (bb-1 ),: )  =  [];  %  deletes  rows  having  svr  assg 
indicator  of  1  (sent  to  svr) 

r  =  size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 
%  here ' s  where  we  implement  h4 : 

IPtemp  =  [IPtemp  zeros (r, 1) ] ;  %  add  a  column  to  hold  the 
current  value 

rows  =  IPtemp (:,1);  %  these  are  the  IP  indices  remaining  in 


IPtemp 

ExpDep  =  ct4(i)  +  es(rows);  %  expected  departure  time  for 
each  IP  if  it  is  sent  to  svc  at  time=ct(i) 

IPtemp ( : , 4 )  = 

(diag (ValueFunction (ExpDep, ValFct Form (rows ) , a (rows ) , par am (rows ,:)))); 
%  computes  h4  for  each  IP  left  in  IPtemp 


IPtemp  =  sortrows ( IPtemp, -4 ) ;  %  sorts  on  h4,  descending 
IPpick  =  IPtemp (1,1);  %choose  IP  w/highest  h4  value 
IPs4 (IPpick, 3)  =  1;  %IP  w/highest  h4  value  goes  to  server 


[v  index]  =  min ( SvrWkld4 ) ;  %  index  is  the  next  avail  svr;  v 
is  the  time  it's  avail. 

SvrAsgn4 ( IPpick, 1 )  =  index;  %assigns  IPpick  to  the  next 
available  server 

if  a (IPpick)  <=  v 

SvrAsgn4 ( IPpick, 2 )  =  v; 

else 

SvrAsgn4 ( IPpick, 2 )  =  a (IPpick);  %  svc  commencement  is 
the  later  of  arrival  time  and  svr  avail  time, 
end 

SvrAsgn4 ( IPpick, 3 )  =  SvrAsgn4 ( IPpick, 2 )  +  s (IPpick);  %  svc 
completion  time 
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server 


is 


SvrWkld4 ( index)  =  SvrAsgn4 (IPpick, 3) ;  %  time  when  this 
next  available. 

%  update  clock  time  (except  for  the  last  iteration) : 
if  i  <  NumPkg 

if  (sum(arr)  >=  i  +  1 )  &&  (min ( SvrWkld4 )  <  a(i+l)) 

ct4(i+l)  =  min ( SvrWkld4 ) ; 

else 


ct4(i+l)  =  max(a(i+l),  min ( SvrWkld4 ) ) ; 

end 


end 


end 

%  save  this  for  output  and  analysis  as  necessary. 

IPdata4  =  [IPs4(:,l:2)  SvrAsgn4  zeros (NumPkg, 1 )] ;  %  col  1  is  IP 
index  (sequenced  by  arrival  time) ; 

%  col  2  is  arrival  time;  col  3  is  server 
%  assigned;  col  4  svc  start  time;  col  5  is  svc 
%  completion  time;  col  6  is  realized  value. 

IPdata4(:,6)  =  diag (ValueFunction (IPdata4 ( : , 5) ,  ValFctForm,  a, 

param) ) ; 


IPdata4sort  =  sortrows ( IPdata4 , 5 ) ;  %  sort  on  svc  completion 

time 

IPdata4sort  =  [IPdata4sort  cumsum ( IPdata4sort ( : , 6 ) ) ] ;  %  add 
additional  col  for  cumulative  realized  value 

cutoff  =  IPdata4sort  (  : , 5)  <=  endtime; 

[g2  h2 ] =min (cutoff ) ; 

RunResults (rep, 2 )  =  IPdata4sort (h2 , 7 ) ;  %  this  is  total  realized 
value  for  h4  for  this  rep 

^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  SERVICE  POLICY  #3:  Highest  Current  Value  (h3) 

9~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

IPs3  =  [(l:NumPkg)'  a  zeros (NumPkg, 1 )] ;  %  first  col  is  IP 
permanent  index 

%  sequenced  by  arrival  time;  second  col  is  arrival  time;  third 
col  is  an 

%  indicator  which  is  1  if  the  IP  has  already  been  sent  to 
server  or  0  if 
%  not . 

SvrWkld3  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time 
assigned  to  each  server. 

SvrAsgn3  =  zeros (NumPkg,  3 ) ;  %each  row  represents  a  package 
(sorted  by  arrival 
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%  time,  just  like  IPs3) ;  first  column 
%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time, 
ct  =  zeros (NumPkg, 1 ) ; 

ct(l)  =  a(l);  %  ct="clock  time"  -  starts  upon  first  arrival. 

%  note  -  ct  is  the  time  an  assignment  decision  is  made  -  at 
%  least  one  server  must  be  available. 

for  i  =  1: NumPkg 

%  extract  from  IPs3  those  pkgs  that  have  1)  arrived  (i.e., 
a<=ct)  and  2) 

%  those  that  have  not  yet  been  sent  to  svr  (3rd  col  of  IPs3 

is  zero) 


IPtemp  =  sortrows (IPs3, 2) ;  %sorts  in  ascending  order  on 

arrival  time 

arr  =  IPtemp (:, 2)  <=  ct(i);  %  arr  =  1  if  IP  has  arrived 
(a<=ct) ;  arr  =  0  otherwise 

%  note:  arr  is  always  size  (NumPkg  x  1) . 

%  Also,  sum(arr)  =  #  of  arrivals  as  of 
%  ct  ( i )  . 

[cc  dd]  =  min (arr);  %  dd  is  the  index  of  the  first  zero 
term,  if  there  is  a  zero 
if  cc  ~=1 

rr=size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp (dd: rr, : )  =  [];  %  deletes  all  rows  w/arrival 
times  greater  than  current  time 
end 

IPtemp  =  sortrows ( IPtemp,  -3);  %  sorts  in  descending  order 

on  col  3 

[aa  bb]  =  min (IPtemp (:, 3) ) ;  %  bb  is  the  row  containing  the 

first  zero 


IPtemp ( 1 : (bb-1 ),: )  =  [];  %  deletes  rows  having  svr  assg 
indicator  of  1  (sent  to  svr) 

r  =  size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp  =  [IPtemp  zeros (r, 1) ] ;  %  add  a  column  to  hold  the 
current  value 

rows  =  IPtemp (:,1);  %  these  are  the  IP  indices  remaining  in 

IPtemp 


IPtemp ( : , 4 )  = 

ValueFunction (ct (i) , ValFctForm (rows) , a (rows) , param (rows, : ) ) ; 

%  computes  the  current  value  of  each  IP  left  in  IPtemp 

IPtemp  =  sortrows ( IPtemp, -4 ) ;  %  sorts  on  current  value, 

descending 

IPpick  =  IPtemp (1,1);  Ichoose  IP  w/highest  current  value 
IPs3 (IPpick, 3)  =  1;  %IP  w/highest  current  value  goes  to 

server 
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[v  index]  =  min ( SvrWkld3 ) ;  %  index  is  the  next  avail  svr;  v 
is  the  time  it's  avail. 

SvrAsgn3 ( IPpick, 1 )  =  index;  %assigns  IPpick  to  the  next 
available  server 

if  a (IPpick)  <=  v 

SvrAsgn3 ( IPpick, 2 )  =  v; 

else 

SvrAsgn3 ( IPpick, 2 )  =  a (IPpick);  %  svc  commencement  is 
the  later  of  arrival  time  and  svr  avail  time, 
end 

SvrAsgn3 (IPpick, 3)  =  SvrAsgn3 ( IPpick, 2 )  +  s (IPpick);  %  svc 
completion  time 

SvrWkld3 ( index)  =  SvrAsgn3 ( IPpick, 3 ) ;  %  time  when  this 
server  is  next  available. 

%  update  clock  time  (except  for  the  last  iteration) : 
if  i  <  NumPkg 

if  (sum(arr)  >=  i  +  1 )  &&  (min ( SvrWkld3 )  <  a(i+l)) 

ct(i+l)  =  min ( SvrWkld3 ) ; 

else 

ct(i+l)  =  max(a(i+l),  min ( SvrWkld3 ) ) ; 

end 


end 


end 

%  save  this  for  output  and  analysis  as  necessary. 

IPdata3  =  [IPs3(:,l:2)  SvrAsgn3  zeros (NumPkg, 1 )] ;  %  col  1  is  IP 
index  (sequenced  by  arrival  time) ; 

%  col  2  is  arrival  time;  col  3  is  server 
%  assigned;  col  4  svc  start  time;  col  5  is  svc 
%  completion  time;  col  6  is  realized  value. 

%  completion  time. 

IPdata3(:,6)  =  diag (ValueFunction (IPdata3 ( : , 5) ,  ValFctForm,  a, 

param) ) ; 


IPdata3sort  =  sortrows ( IPdata3 , 5 ) ;  %  sort  on  svc  completion 

time 

IPdata3sort  =  [IPdata3sort  cumsum ( IPdata3sort ( : , 6) ) ] ;  %  add 
additional  col  for  cumulative  realized  value 

cutoff  =  IPdata3sort ( : , 5)  <=  endtime; 

[g2  h2 ] =min (cutoff ) ; 

RunResults (rep, 3)  =  IPdata3sort (h2 , 7 ) ;  %  this  is  total  realized 
value  for  h3  for  this  rep 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  SERVICE  POLICY  #7:  LIFO  -  h7 

9~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 
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first  col  is  IP 


IPs7  =  [(l:NumPkg)'  a  zeros (NumPkg, 1 )] ;  % 
permanent  index 

%  sequenced  by  arrival  time;  second  col  is  arrival  time;  third 
col  is  an 

%  indicator  which  is  1  if  the  IP  has  already  been  sent  to 
server  or  0  if 
%  not . 

SvrWkld7  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time 
assigned  to  each  server. 

SvrAsgn7  =  zeros (NumPkg, 3 ) ;  %each  row  represents  a  package 
(sorted  by  arrival 

%  time,  just  like  IPs);  first  column 
%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time. 
ct7  =  zeros (NumPkg, 1 ) ; 

ct7(l)  =  a(l);  %  ct="clock  time"  -  starts  upon  first  arrival. 

%  note  -  ct  is  the  time  an  assignment  decision  is  made  -  at 
%  least  one  server  must  be  available. 

for  i  =  1: NumPkg 

%  extract  from  IPs  those  pkgs  that  have  1)  arrived  (i.e., 
a<=ct)  and  2) 

%  those  that  have  not  yet  been  sent  to  svr  (3rd  col  of  IPs 

is  zero) 


IPtemp  =  sortrows ( IPs7 , 2 ) ;  %sorts  in  ascending  order  on 

arrival  time 

arr  =  IPtemp (:, 2)  <=  ct7(i);  %  arr  =  1  if  IP  has  arrived 
(a<=ct7);  arr  =  0  otherwise 

%  note:  arr  is  always  size  (NumPkg  x  1) . 

%  Also,  sum(arr)  =  #  of  arrivals  as  of 
%  ct  ( i )  . 

[cc  dd]  =  min (arr);  %  dd  is  the  index  of  the  first  zero 
term,  if  there  is  a  zero 

if  cc  ~=1  %  ie,  if  there  is  at  least  one  IP  that  hasn't 

"arrived" 

rr=size (IPtemp, 1) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp (dd: rr, : )  =  [];  %  deletes  all  rows  w/arrival 
times  greater  than  current  time 
end 

IPtemp  =  sortrows ( IPtemp,  -3);  %  sorts  in  descending  order 
on  col  3  (sent  to  svr) 

[aa  bb]  =  min (IPtemp (:, 3) ) ;  %  bb  is  the  row  containing  the 

first  zero 


IPtemp (1 : (bb-1) ,: )  =  [];  %  deletes  rows  having  svr  assg 
indicator  of  1  (sent  to  svr) 

r  =  size ( IPtemp,  1 ) ;  %  number  of  rows  remaining  in  IPtemp 
%  here ' s  where  we  implement  h7 : 
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[ar  lastin]  =  max ( IPtemp ( : , 2 ) ) ; 

IPpick  =  IPtemp (lastin, 1) ; 

IPs7 (IPpick, 3)  =  1;  %IP  w/highest  h7  value  goes  to  server 

[v  index]  =  min ( SvrWkld7 ) ;  %  index  is  the  next  avail  svr;  v 
is  the  time  it's  avail. 

SvrAsgn7 ( IPpick, 1 )  =  index;  %assigns  IPpick  to  the  next 
available  server 

if  a (IPpick)  <=  v 

SvrAsgn7 ( IPpick, 2 )  =  v; 

else 

SvrAsgn7 ( IPpick, 2 )  =  a (IPpick);  %  svc  commencement  is 
the  later  of  arrival  time  and  svr  avail  time. 

end 

SvrAsgn7 ( IPpick, 3 )  =  SvrAsgn7 ( IPpick, 2 )  +  s (IPpick);  %  svc 
completion  time 

SvrWkld7 ( index)  =  SvrAsgn7 ( IPpick, 3 ) ;  %  time  when  this 
server  is  next  available. 


%  update  clock  time  (except  for  the  last  iteration) : 
if  i  <  NumPkg 

if  (sum(arr)  >=  i  +  1 )  &&  (min ( SvrWkld7 )  <  a(i+l)) 

ct7(i+l)  =  min ( SvrWkld7 ) ; 

else 


ct7(i+l)  =  max(a(i+l),  min ( SvrWkld7 ) ) ; 

end 


end 


end 

%  save  this  for  output  and  analysis  as  necessary. 

IPdata7  =  [IPs7(:,l:2)  SvrAsgn7  zeros (NumPkg, 1 )] ;  %  col  1  is  IP 
index  (sequenced  by  arrival  time) ; 

%  col  2  is  arrival  time;  col  3  is  server 
%  assigned;  col  4  svc  start  time;  col  5  is  svc 
%  completion  time;  col  6  is  realized  value. 

IPdata7(:,6)  =  diag (ValueFunction (IPdata7 ( : , 5) ,  ValFctForm,  a, 

param) ) ; 


IPdata7sort  =  sortrows ( IPdata7 , 5 ) ;  %  sort  on  svc  completion 

time 

IPdata7sort  =  [IPdata7sort  cumsum ( IPdata7sort ( : , 6) ) ] ;  %  add 
additional  col  for  cumulative  realized  value 

cutoff  =  IPdata7sort ( : , 5)  <=  endtime; 

[g2  h2 ] =min (cutoff ) ; 

RunResults (rep, 4 )  =  IPdata7sort (h2, 7) ;  %  this  is  total  realized 
value  for  LIFO  for  this  rep 
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9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  SERVICE  POLICY  #8:  Random  -  h8 

9~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

IPs8  =  [(l:NumPkg)'  a  zeros (NumPkg, 1 )] ;  %  first  col  is  IP 
permanent  index 

%  sequenced  by  arrival  time;  second  col  is  arrival  time;  third 
col  is  an 

%  indicator  which  is  1  if  the  IP  has  already  been  sent  to 
server  or  0  if 
%  not . 

SvrWkld8  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time 
assigned  to  each  server. 

SvrAsgn8  =  zeros (NumPkg, 3 ) ;  %each  row  represents  a  package 
(sorted  by  arrival 

%  time,  just  like  IPs);  first  column 
%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time. 
ct8  =  zeros (NumPkg, 1 ) ; 

ct8(l)  =  a(l);  %  ct="clock  time"  -  starts  upon  first  arrival. 

%  note  -  ct  is  the  time  an  assignment  decision  is  made  -  at 
%  least  one  server  must  be  available. 

for  i  =  1: NumPkg 

%  extract  from  IPs  those  pkgs  that  have  1)  arrived  (i.e., 
a<=ct)  and  2) 

%  those  that  have  not  yet  been  sent  to  svr  (3rd  col  of  IPs 

is  zero) 


IPtemp  =  sortrows ( IPs8 , 2 ) ;  %sorts  in  ascending  order  on 

arrival  time 

arr  =  IPtemp (:, 2)  <=  ct8(i);  %  arr  =  1  if  IP  has  arrived 
(a<=ct8) ;  arr  =  0  otherwise 

%  note:  arr  is  always  size  (NumPkg  x  1) . 

%  Also,  sum(arr)  =  #  of  arrivals  as  of 
%  ct  ( i )  . 

[cc  dd]  =  min (arr);  %  dd  is  the  index  of  the  first  zero 
term,  if  there  is  a  zero 

if  cc  ~=1  %  ie,  if  there  is  at  least  one  IP  that  hasn't 

"arrived" 

rr=size (IPtemp, 1) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp (dd: rr, : )  =  [];  %  deletes  all  rows  w/arrival 
times  greater  than  current  time 
end 

IPtemp  =  sortrows ( IPtemp,  -3);  %  sorts  in  descending  order 
on  col  3  (sent  to  svr) 

[aa  bb]  =  min (IPtemp (:, 3) ) ;  %  bb  is  the  row  containing  the 

first  zero 


IPtemp (1 : (bb-1) ,: )  =  [];  %  deletes  rows  having  svr  assg 
indicator  of  1  (sent  to  svr) 
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r  =  size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 

%  here's  where  we  implement  h8  -  pick  an  IP  from  IPtemp  at 
random  to  send  to  svr: 

pick  =  ceil ( r*rand) ; 

IPpick  =  IPtemp (pick, 1) ; 

IPs8 (IPpick, 3)  =  1;  %IP  w/highest  h8  value  goes  to  server 

[v  index]  =  min ( SvrWkld8 ) ;  %  index  is  the  next  avail  svr;  v 
is  the  time  it's  avail. 

SvrAsgn8 ( IPpick, 1 )  =  index;  %assigns  IPpick  to  the  next 
available  server 

if  a (IPpick)  <=  v 

SvrAsgn8 ( IPpick, 2 )  =  v; 

else 

SvrAsgn8 ( IPpick, 2 )  =  a (IPpick);  %  svc  commencement  is 
the  later  of  arrival  time  and  svr  avail  time. 

end 

SvrAsgn8 (IPpick, 3)  =  SvrAsgn8 ( IPpick, 2 )  +  s (IPpick);  %  svc 
completion  time 

SvrWkld8 ( index)  =  SvrAsgn8 (IPpick, 3) ;  %  time  when  this 
server  is  next  available. 


%  update  clock  time  (except  for  the  last  iteration) : 
if  i  <  NumPkg 


if  (sum(arr) 
ct8 ( i+1 ) 

else 

ct8 ( i+1 ) 

end 


>=  i  +  1 )  &&  (min ( SvrWkld8 )  <  a (i+1)) 

=  min ( SvrWkld8 ) ; 

=  max (a (i+1),  min ( SvrWkld8 ) ) ; 


end 


end 

%  save  this  for  output  and  analysis  as  necessary. 

IPdata8  =  [IPs8(:,l:2)  SvrAsgn8  zeros (NumPkg, 1 )] ;  %  col  1  is  IP 
index  (sequenced  by  arrival  time) ; 

%  col  2  is  arrival  time;  col  3  is  server 
%  assigned;  col  4  svc  start  time;  col  5  is  svc 
%  completion  time;  col  6  is  realized  value. 

IPdata8(:,6)  =  diag (ValueFunction (IPdata8 ( : , 5) ,  ValFctForm,  a, 

param) ) ; 


IPdata8sort  =  sortrows ( IPdata8 , 5 ) ;  %  sort  on  svc  completion 

time 

IPdata8sort  =  [IPdata8sort  cumsum ( IPdata8sort ( : , 6) ) ] ;  %  add 
additional  col  for  cumulative  realized  value 

cutoff  =  IPdata8sort ( : , 5)  <=  endtime; 
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[g2  h2 ] =min (cutoff ) ; 


RunResults (rep, 5)  =  IPdata8sort (h2 , 7 ) ;  %  this  is  total  realized 
value  for  RANDOM  for  this  rep 

^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  SERVICE  POLICY  #9:  "Status  Quo" 

9~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


%  this  heuristic  attempts  to  replicate  the  current  strategy  for 

choosing 

%  IPs  to  process  by  combining  three  of  the  above  heuristics: 
h3 ,  LIFO, 

%  and  Random.  These  can  be  "evenly"  combined,  or  weighted, 
weight  =  [.3  .6  .1];  %  weights  for  h3,  LIFO,  and  Random, 
respectively . 

%Must  sum  to  one. 
st  =  rand (NumPkg, 1 ) ; 
strategy  =  zeros (NumPkg, 1 ) ; 
str  =  [0  cumsum (weight) ] ; 
for  j  =1 : 3 

strategy  =  strategy  +  j*(st  >  str(j)  &  st  <=  str(j+l)); 

end 

%  strategy  is  a  vector  of  length  NumPkg  that  indicates  a  1  if 
h3  is  to  be 

%  used,  2  if  LIFO  is  to  be  used,  and  3  if  Random  is  to  be  used. 

Its 


%  elements  are  randomly  generated  according  to  the  vector 
"weight . " 


IPs9  =  [(l:NumPkg)'  a  zeros (NumPkg, 1 )] ;  %  first  col  is  IP 
permanent  index 

%  sequenced  by  arrival  time;  second  col  is  arrival  time;  third 
col  is  an 

%  indicator  which  is  1  if  the  IP  has  already  been  sent  to 
server  or  0  if 
%  not . 

SvrWkld9  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time 
assigned  to  each  server. 

SvrAsgn9  =  zeros (NumPkg, 3 ) ;  %each  row  represents  a  package 
(sorted  by  arrival 

%  time,  just  like  IPs);  first  column 
%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time. 
ct9  =  zeros (NumPkg, 1 ) ; 

ct9(l)  =  a(l);  %  ct="clock  time"  -  starts  upon  first  arrival. 

%  note  -  ct  is  the  time  an  assignment  decision  is  made  -  at 
%  least  one  server  must  be  available. 

for  i  =  1: NumPkg 

%  extract  from  IPs  those  pkgs  that  have  1)  arrived  (i.e., 
a<=ct)  and  2) 
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is  zero) 


those  that  have  not  yet  been  sent  to  svr  (3rd  col  of  IPs 


IPtemp  =  sortrows ( IPs9, 2 ) ;  %sorts  in  ascending  order  on 

arrival  time 

arr  =  IPtemp (:, 2)  <=  ct9(i);  %  arr  =  1  if  IP  has  arrived 
(a<=ct9) ;  arr  =  0  otherwise 

%  note:  arr  is  always  size  (NumPkg  x  1) . 

%  Also,  sum(arr)  =  #  of  arrivals  as  of 
%  ct  ( i )  . 

[cc  dd]  =  min (arr);  %  dd  is  the  index  of  the  first  zero 
term,  if  there  is  a  zero 

if  cc  ~=1  %  ie,  if  there  is  at  least  one  IP  that  hasn't 

"arrived" 

rr=size (IPtemp, 1) ;  %  number  of  rows  remaining  in  IPtemp 
IPtemp (dd: rr, : )  =  [];  %  deletes  all  rows  w/arrival 
times  greater  than  current  time 
end 

IPtemp  =  sortrows ( IPtemp,  -3);  %  sorts  in  descending  order 
on  col  3  (sent  to  svr) 

[aa  bb]  =  min (IPtemp (:, 3) ) ;  %  bb  is  the  row  containing  the 

first  zero 


IPtemp (1 : (bb-1) ,: )  =  [];  %  deletes  rows  having  svr  assg 
indicator  of  1  (sent  to  svr) 

r  =  size ( IPtemp, 1 ) ;  %  number  of  rows  remaining  in  IPtemp 

%  here's  where  we  implement  h9:  based  on  value  of 
strategy (i) ,  choose 

%  h3,  LIFO,  or  Random. 

if  strategy (i)  ==  3  %  random 
pick  =  ceil (r*rand) ; 

IPpick  =  IPtemp (pick, 1) ; 
elseif  strategy(i)  ==  2  %  LIFO 

[ar  lastin]  =  max ( IPtemp (:, 2 )) ; 

IPpick  =  IPtemp (lastin, 1) ; 
else  %  h3 

IPtemp  =  [IPtemp  zeros (r,  1) ] ;  %  add  a  column  to  hold 
the  current  value 

rows  =  IPtemp (:,1);  %  these  are  the  IP  indices 
remaining  in  IPtemp 


IPtemp ( : , 4 )  = 

ValueFunction (ct ( i ) , ValFctForm ( rows ) , a (rows ) , par am (rows , : ) ) ; 

%  computes  the  current  value  of  each  IP  left  in  IPtemp 


descending 

value 


IPtemp  =  sortrows ( IPtemp, -4 ) ;  %  sorts  on  current  value, 
IPpick  =  IPtemp (1,1);  %choose  IP  w/highest  current 


end 
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IPs9 (IPpick, 3)  =  1;  %IP  w/highest  h9  value  goes  to  server 

[v  index]  =  min ( SvrWkld9 ) ;  %  index  is  the  next  avail  svr;  v 
is  the  time  it's  avail. 

SvrAsgn9 ( IPpick, 1 )  =  index;  %assigns  IPpick  to  the  next 
available  server 

if  a (IPpick)  <=  v 

SvrAsgn9 ( IPpick, 2 )  =  v; 

else 

SvrAsgn9 ( IPpick, 2 )  =  a (IPpick);  %  svc  commencement  is 
the  later  of  arrival  time  and  svr  avail  time. 

end 

SvrAsgn9 (IPpick, 3)  =  SvrAsgn9 ( IPpick, 2 )  +  s (IPpick);  %  svc 

completion  time 

SvrWkld9 ( index)  =  SvrAsgn9 (IPpick, 3) ;  %  time  when  this 
server  is  next  available. 


%  update  clock  time  (except  for  the  last  iteration) : 
if  i  <  NumPkg 

if  (sum(arr)  >=  i  +  1 )  &&  (min (SvrWkld9)  <  a(i+l)) 

ct9(i+l)  =  min (SvrWkld9) ; 

else 


ct9(i+l)  =  max(a(i+l),  min (SvrWkld9) ) ; 

end 


end 


end 

%  save  this  for  output  and  analysis  as  necessary. 

IPdata9  =  [IPs9(:,l:2)  SvrAsgn9  zeros (NumPkg, 1 )] ;  %  col  1  is  IP 
index  (sequenced  by  arrival  time) ; 

%  col  2  is  arrival  time;  col  3  is  server 
%  assigned;  col  4  svc  start  time;  col  5  is  svc 
%  completion  time;  col  6  is  realized  value. 

IPdata9(:,6)  =  diag (ValueFunction (IPdata9 ( : , 5) ,  ValFctForm,  a, 

param) ) ; 


IPdata9sort  =  sortrows (IPdata9, 5) ;  %  sort  on  svc  completion 

time 

IPdata9sort  =  [IPdata9sort  cumsum ( IPdata9sort ( : , 6 ) ) ] ;  %  add 
additional  col  for  cumulative  realized  value 

cutoff  =  IPdata9sort ( : , 5)  <=  endtime; 

[g2  h2 ] =min (cutoff ) ; 

RunResults (rep, 6)  =  IPdata9sort (h2 , 7 ) ;  %  this  is  total  realized 
value  for  StatusQuo  for  this  rep 

^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

%  SERVICE  POLICY  #1:  FIFO: 
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9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


SvrWkld  =  zeros (NumSvr, 1) ;  %  tracks  total  service  time  assigned 
to  each  server. 

SvrAsgn  =  zeros (NumPkg, 3 ) ;  %each  row  represents  a  package; 
first  column 

%is  the  server  assigned;  second  column  is 
%processing  start  time;  third  column  is 
%processing  completion  time. 

%NOTE :  the  following  construction  assumes  that  pkgs  are  listed 

in  order  of 

%arrival  time! ! 
for  j  =1 : NumPkg 

[v  index]  =  min ( SvrWkld) ;  %  index  gives  the  next  server 

available 

SvrAsgn (j,l)  =  index;  %assign  the  next  available  server  to 

pkg  j  . 

if  a ( j )  <=  v 

SvrAsgn (j, 2)  =  v; 

else 

SvrAsgn (j, 2)  =  a(j);  %  service  commencement  is  the 
smallest  of  v  and  a(j) 
end 

SvrAsgn (j, 3)  =  SvrAsgn (j, 2)  +  s(j);  %  computes  service 
completion  time  for  pkg  j . 

SvrWkld ( index)  =  SvrAsgn ( j , 3 ) ;  %  time  when  this  server  is 
next  available 

end 

val  r  =  diag (ValueFunction ( SvrAsgn ( : , 3 ) , ValFctForm, a, param) ) ; 
%realized  value  for  each  package 

%  the  sorting  below  accounts  for  the  fact  that  even  though  pkgs 

are 

%  ENTERING  service  in  order  of  arrival,  they  do  not  necessarily 

exit  the 

%  server  in  the  same  order  (some  later  arriving  pkgs  are 

finished 

%  processing  before  some  earlier  arriving  ones) . 
valmat  =  [ SvrAsgn (:, 3 )  val  r] ;  %  svc  completion  times  paired 
w/realized  values 

valmatsort  =  sortrows (valmat, 1) ;  %  sort  on  svc  completion  times 

val  r  t  =  cumsum (valmatsort (:, 2 )) ;  %  cumulatively  sum  the 
realized  values  for  each  pkg 

SvrCompleteTimes  =  valmatsort (:, 1 ) ;  %NOTE :  SvrCompleteTimes  is 
NOT  sorted  in  pkg  number  order! ! 

%plot (SvrCompleteTimes, val  r  t)  %  plots  realized  value  vs.  time 
for  all  pkgs 
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cut of f=SvrCompleteTimes<=endtime; 

[g  h] =min (cutoff ) ;  %  h  is  the  index  of  the  first  zero,  which 

equates  to  the  first  svc  complete  time  later  than  endtime. 


RunResults (rep, 7 )  =  val  r  t (h) ;  %  this  is  total  realized  value 
for  FIFO  for  this  rep 

end 

meanrow=NumReps  +  l  ; 
stdrow=NumReps+2  ; 

TempResults  =  RunResults ( 1 : NumReps ,:)  ; 

RunResults (meanrow, : )  =  Mean (TempResults) ;  %mean  for  each  policy 
RunResults ( stdrow, : )  =  std (TempResults) ;  %  standard  deviation 


overall  mean (run,  1)  =  parameters (run, 1) ;  %  run  number 
overall_std (run, 1 )  =  parameters (run, 1) ; 
overall  mean (run, 2 : 8)  =  RunResults (meanrow, :) ; 
overall  std (run, 2: 8)  =  RunResults ( stdrow, :) ; 


csvwrite ( [ ' run  ' , num2str (parameters ( run, 1)),'  results . csv '] , RunResults ) 


end 

toe 

csvwrite ( [num2str (parameters (1, 1) ), '  overall  mean . csv ']? overall  mean); 
csvwrite ( [num2str (parameters (1,1) ) , ' _overall_std . csv ' ] , overall_std) ; 
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E. 


DYNAMIC  SIMULATION  RUN  MATRIX 


ValFctMix  Alpha  Beta  Beta 
(E/L/S)  (lin.exp)  (lin.exp)  (step) 
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75/15/10 


75/15/10 
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75/15/10 


75/15/10 


75/15/10 
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75/15/10 
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75/15/10 
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0.631, 0.5,  0.8 
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[0.631, 0.5,  0.8] 


1.2,  0.2,  0.6 
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0.631,  0.5,  0.8 
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[0.631, 0.5,  0.8] 
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[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 
[0.134,  0.1,  1] 


Run 

Reps 

#  IPs 

Num- 

Class 

Num- 

Svrs 

ValFctMix 

(E/L/S) 

Alpha 

(lin.exp) 

Beta 

(lin.exp) 

Beta 

(step) 

tau 

(step) 

Lambda 

Mu 

51 

250 

1000 

3 

3 

75/15/10 

H 

L 

L 

L 

(0. 631, 0.5,  0.81 

[0.134,  0.1, 

11 

52 

250 

1000 

3 

3 

75/15/10 

H 

L 

L 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

1] 

53 

250 

1000 

3 

3 

75/15/10 

H 

L 

H 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

54 

250 

1000 

3 

3 

75/15/10 

H 

L 

H 

H 

fl  .2,  0.2,  0.61 

[0.134,  0.1, 

11 

55 

250 

1000 

3 

3 

75/15/10 

H 

L 

H 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

56 

250 

1000 

3 

3 

75/15/10 

H 

L 

H 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

57 

250 

1000 

3 

3 

75/15/10 

H 

H 

L 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

58 

250 

1000 

3 

3 

75/15/10 

H 

H 

L 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

1] 

59 

250 

1000 

3 

3 

75/15/10 

H 

H 

L 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

60 

250 

1000 

3 

3 

75/15/10 

H 

H 

L 

H 

fl  .2,  0.2,  0.61 

[0.134,  0.1, 

11 

61 

250 

1000 

3 

3 

75/15/10 

H 

H 

H 

L 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

1] 

62 

250 

1000 

3 

3 

75/15/10 

H 

H 

H 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

63 

250 

1000 

3 

3 

75/15/10 

H 

H 

H 

H 

[0.631,  0.5,  0.81 

[0.134,  0.1, 

11 

64 

250 

1000 

3 

3 

75/15/10 

H 

H 

H 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

1] 

65 

250 

1000 

3 

5 

33/33/33 

L 

L 

L 

L 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

66 

250 

1000 

3 

5 

33/33/33 

L 

L 

L 

L 

]1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

67 

250 

1000 

3 

5 

33/33/33 

L 

L 

L 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

1] 

68 

250 

1000 

3 

5 

33/33/33 

L 

L 

L 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

69 

250 

1000 

3 

5 

33/33/33 

L 

L 

H 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

70 

250 

1000 

3 

5 

33/33/33 

L 

L 

H 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

1] 

71 

250 

1000 

3 

5 

33/33/33 

L 

L 

H 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

72 

250 

1000 

3 

5 

33/33/33 

L 

L 

H 

H 

]1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

73 

250 

1000 

3 

5 

33/33/33 

L 

H 

L 

L 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

1] 

74 

250 

1000 

3 

5 

33/33/33 

L 

H 

L 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

75 

250 

1000 

3 

5 

33/33/33 

L 

H 

L 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

76 

250 

1000 

3 

5 

33/33/33 

L 

H 

L 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

1] 

77 

250 

1000 

3 

5 

33/33/33 

L 

H 

H 

L 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

78 

250 

1000 

3 

5 

33/33/33 

L 

H 

H 

L 

]1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

79 

250 

1000 

3 

5 

33/33/33 

L 

H 

H 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

1] 

80 

250 

1000 

3 

5 

33/33/33 

L 

H 

H 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

81 

250 

1000 

3 

5 

33/33/33 

H 

L 

L 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

82 

250 

1000 

3 

5 

33/33/33 

H 

L 

L 

L 

[1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

83 

250 

1000 

3 

5 

33/33/33 

H 

L 

L 

H 

[0.631, 0.5,  0.8] 

o 

CO 

o 

11 

84 

250 

1000 

3 

5 

33/33/33 

H 

L 

L 

H 

]1.2,  0.2,  0.61 

o 

CO 

CD 

11 

85 

250 

1000 

3 

5 

33/33/33 

H 

L 

H 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

86 

250 

1000 

3 

5 

33/33/33 

H 

L 

H 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

87 

250 

1000 

3 

5 

33/33/33 

H 

L 

H 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

88 

250 

1000 

3 

5 

33/33/33 

H 

L 

H 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

89 

250 

1000 

3 

5 

33/33/33 

H 

H 

L 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

90 

250 

1000 

3 

5 

33/33/33 

H 

H 

L 

H 

]1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

91 

250 

1000 

3 

5 

33/33/33 

H 

H 

L 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

92 

250 

1000 

3 

5 

33/33/33 

H 

H 

L 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

93 

250 

1000 

3 

5 

33/33/33 

H 

H 

H 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

94 

250 

1000 

3 

5 

33/33/33 

H 

H 

H 

H 

[1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

95 

250 

1000 

3 

5 

33/33/33 

H 

H 

H 

L 

[0.631, 0.5,  0.8] 

[0.134,  0.1, 

11 

96 

250 

1000 

3 

5 

33/33/33 

H 

H 

H 

L 

]1.2,  0.2,  0.61 

[0.134,  0.1, 

11 

97 

250 

1000 

3 

5 

75/15/10 

L 

L 

L 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

98 

250 

1000 

3 

5 

75/15/10 

L 

L 

L 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

11 

99 

250 

1000 

3 

5 

75/15/10 

L 

L 

L 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1, 

11 

100 

250 

1000 

3 

5 

75/15/10 

L 

L 

L 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1, 

1] 

134 


Run 

Reps 

#  IPs 

Num- 

Num- 

ValFctMix 

Alpha 

Beta 

Beta 

tau 

Lambda 

Mu 

Class 

Svrs 

(E/L/S) 

(lin.exp) 

(lin.exp) 

(step) 

(step) 

inn 

■IIMf 

3 

5 

75/15/10 

L 

L 

H 

H 

[0.631, 0.5,  0.8] 

[0.134,  0.1,  1] 

102 

250 

1000 

3 

5 

75/15/10 

L 

L 

H 

H 

[1.2,  0.2,  0.61 

mm 

103 

250 

1000 

3 

5 

75/15/10 

L 

L 

H 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

104 

1000 

3 

5 

75/15/10 

L 

L 

H 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1,  1] 

105 

1000 

3 

5 

75/15/10 

L 

H 

L 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

106 

1000 

3 

5 

75/15/10 

L 

H 

L 

H 

[1.2,  0.2,  0.61 

Rmn 

107 

1  250  | 

1000 

3 

5 

75/15/10 

L 

H 

L 

L 

[0.631, 0.5,  0.8] 

|  [0.134,  0.1,  11  | 

108 

1000 

3 

5 

75/15/10 

L 

H 

L 

L 

[1.2,  0.2,  0.61 

Rmn 

109 

1  250  | 

1000 

3 

5 

75/15/10 

L 

H 

H 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

110 

1000 

3 

5 

75/15/10 

L 

H 

H 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1,  1] 

111 

1000 

3 

5 

75/15/10 

L 

H 

H 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

112 

1000 

3 

5 

75/15/10 

L 

H 

H 

H 

[1.2,  0.2,  0.61 

Rmn 

i/TEl 

3 

5 

75/15/10 

H 

L 

L 

L 

[0.631, 0.5,  0.8] 

|  [0.134,  0.1,  11  | 

114 

250 

1000 

3 

5 

75/15/10 

H 

L 

L 

L 

[1.2,  0.2,  0.61 

Rmn 

115 

250 

1000 

3 

5 

75/15/10 

H 

L 

L 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

116 

1000 

3 

5 

75/15/10 

H 

L 

L 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1,  1] 

117 

1000 

3 

5 

75/15/10 

H 

L 

H 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

118 

1000 

3 

5 

75/15/10 

H 

L 

H 

L 

[1.2,  0.2,  0.61 

Rmn 

119 

1  250  | 

1000 

3 

5 

75/15/10 

H 

L 

H 

H 

[0.631, 0.5,  0.8] 

|  [0.134,  0.1,  11  | 

120 

1000 

3 

5 

75/15/10 

H 

L 

H 

H 

[1.2,  0.2,  0.61 

Rmn 

121 

1  250  | 

1000 

3 

5 

75/15/10 

H 

H 

L 

L 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

3 

5 

75/15/10 

H 

H 

L 

L 

[1.2,  0.2,  0.6] 

[0.134,  0.1,  1] 

123 

1000 

3 

5 

75/15/10 

H 

H 

L 

H 

[0.631, 0.5,  0.81 

[0.134,  0.1,  11 

124 

1000 

3 

5 

75/15/10 

H 

H 

L 

H 

[1.2,  0.2,  0.61 

Rmn 

1E3 

■HMf 

3 

5 

75/15/10 

H 

H 

H 

L 

[0.631, 0.5,  0.8] 

|  [0.134,  0.1,  11  | 

126 

1  250  | 

1000 

3 

5 

75/15/10 

H 

H 

H 

L 

[1.2,  0.2,  0.61 

Rmn 

127 

E5I 

1000 

3 

5 

75/15/10 

H 

H 

H 

H 

[0.631,  0.5,  0.81 

[0.134,  0.1,  11 

128 

1  250  | 

1000 

3 

5 

75/15/10 

H 

H 

H 

H 

[1.2,  0.2,  0.6] 

[0.134,  0.1,  1] 

135 
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APPENDIX  B.  DISCRETE  STOCHASTIC  INFORMATION  FLOW 

MODEL  MATLAB  CODE 


%%  This  m-file  simulates  the  system  x(k+l)  =  Ax(k)  +  Bu(k),  where  u(k) 
is 

%%  a  vector  of  random  variables  depicting  1)  the  number  of  info  pkgs 
%%  arriving  in  Oval  2  during  time  period  k,  and  2)  which  processing 
time 

%%  "bin"  these  pkgs  are  assigned  to. 

%%  Version  1.0,  24  Jan  2007 
%%  Donovan  Phillips 

clear  all 

%  First,  define  some  parameters: 
n  =  1000;  %  number  of  runs  in  the  simulation 

p  =  51;  %  number  of  data  pkgs  to  generate  each  run  -  same  number  as 

orig.  data 

di  =  2;  %  delta  i:  size  (in  min)  of  each  processing  time  "bin" 

dt  =  2 ;  %  time  step  size  (min) 

mu  arrtime  =  0.843;  %  average  pkg  inter-arrival  time 

mu  svctime  =  7.03529;  %  average  pkg  processing  time 

max  endtime  =  300;  %  simulation  end  time,  min 
alpha  svc  =  0.410298; 

beta  svc  =  17.1468;  %  parameters  for  beta  fit  of  svc  time  data 
alpha  arr  =  0.7351; 

beta  arr  =  1.614;  %  parameters  for  beta  fit  of  arr  time  data 

%%  initialize  the  output  matrix: 
vol  =  zeros (n, max  endtime/dt); 
t  =  dt:dt:max  endtime; 

%  Now,  perform  the  simulation: 

for  i=l : n 

rand (' state ', sum ( 100*clock) ) ;  %randomly  sets  the  random  number 
generator  state  (or  "seed") 

%inter  arr  times  =  exprnd (mu  arrtime, p, 1) ; 

inter  arr  times  =  gamrnd(alpha  arr,  beta  arr,  p,l); 

rand ( ' state ' , sum (100*clock)  )  ; 

%svc  times  =  exprnd (mu  svctime, p, 1) ; 

svc_times  =  gamrnd(alpha  svc,  beta  svc,  p,  1) ; 

arr_times  =  zeros (p,l); 

arr  times (1)  =  inter  arr  times (1); 

for  j=l : (p-1) 

arr  times (j+1)  =  arr  times (j)  +  inter  arr  times (j+1); 

end 

P=zeros (p, 3 ) ; 
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P ( : , 1 ) =arr_times ; 

P ( : , 3) =svc_times; 

P(:,2)=arr  times  +  svc  times; 

d  =  CoefMatrixMod (di, dt, P) ;  %  rows  of  d  are  processing  time  bins; 
cols  are  time  steps 
s  =  size (d) ; 

si  =  s(l);  %  number  of  rows  in  d 

s2  =  s (2) ;  %  number  of  cols  in  d;  note:  size  of  d  will  vary  with 
each  run. 

bins  =  (di : di :  (sl*2 ) )  '  ; 
for  k=l:s2 

vol(i,k)  =  sum (bins . *d (:, k) ) ;  %computes  the  "volume"  of  info  in 
oval  2  at  time  period  k 
end 

end 

sumvol  =  sum(vol); 
close  all 
plot (t, vol) 
hold  on 

ad  =  importdata ( ' processtimeNoIntelDump22 . csv ' ) ;  %actual  data,  for 
comparison 

matrix  =  CoefMatrixMod (di, dt, ad) ; 
bins  =  (di:di:(2*size (matrix,  1) ) )  '  ; 

vol  actual  =  zeros (1, max  endtime/dt); 

for  i=l : (size (matrix, 2) ) 

vol_actual ( i ) =sum (bins . *matrix ( : ,  i ) ) ; 

end 

plot ( t , vol_actual , ' kx : ' ) 
xlabel ( ' Time  (min) ' ) ; 

ylabel (' Total  Volume  in  Oval  2  (min) '); 

xlim([0  1 . 05*dt*max (max ( find ( sumvol ))  ,  max (find (vol  actual)))]); 

hold  off 
figure 
hold  on 

avg  vol  =  mean (vol); 

stdev_vol  =  std(vol); 

ci  =  1 . 0*stdev_vol; 

upper  =  avg  vol  +  ci; 

lower  =  max (avg  vol  -  ci,0); 

plot (t, avg_vol) 

plot (t, upper, 'b--') 
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plot (t, lower, 'b--') 
plot ( t , vol_actual , ' kx : ' ) 
xlabel ( ' Time  (min) ' ) ; 

ylabel (' Total  Volume  in  Oval  2  (min) '); 

xlim([0  1 . 05*dt*max (max (find (sumvol) ) ,  max(find(vol  actual)))]) 
hold  off 
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